Notebooks¶
The following tutorials and research applications refer to Jupyter notebooks provided in the repository. These notebooks can be started directly with Binder
Tutorial 1: Introduction¶
In this tutorial, the basic steps of using pyLabFEA for elastic materials are demonstrated.
Finite Element Tutorial in Python¶
Introduction to the Python Laboratory for Finite Element Analysis (pyLabFEA). This tutorial covers FEA with elements in one or two dimensions to evaluate the mechanical behavior of laminate geometries with linear elastic materials. See the online documentation of the pyLabFEA package for detailed information. NumPy (http://www.numpy.org) is used for mathematical operations on arrays and matplotlib (https://matplotlib.org/) for the visualization of numerical results.
Author: Alexander Hartmaier, ICAMS / Ruhr-Universität Bochum
March 2020
This work is licensed under a Creative
Commons Attribution-NonCommercial-ShareAlike 4.0 International License
(CC-BY-NC-SA)
The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)
Python Laboratory for Finite Element Analysis (pyLabFEA)¶
In the first step, the pyLabFEA package, containing the modules model and material, is imported.
import pylabfea.model as FE
from pylabfea.material import Material
Model¶
The class Model is invoked to create an objecte for the finite element model to be used. This object contains all attributes of the model and defines the methods. To use the model for FEA, the following steps are performed: (i) geometry definition, (ii) material definition, (iii) applying boundary conditions, (iv) meshing, (v) solving of linear equations, and (vi) calculation of element stresses and strains by calling the corresponding methods of the class Model.
Parameters:
dim (optional): dim=1 for 1-d model, dim=2 for 2-d model (default: dim=1)
planestress (optional, used only in 2-d models): True if plane stress conditions shall be considered (default: planestress=False)
Methods:
geom: define (laminate) geometry with sections
assign: assign material properties to sections
bcleft, bcright, bctop, bcbot: define boundary conditions
mesh: generate mesh (1-d elements with linear or quadratic shape functions or 2-d mesh with linear quadriliteratl elements)
solve: setup stiffness matrix and calculate displacements satisfying mechanical equilibrium, i.e. $f=0$ for internal nodes, for given boundary conditions
calc_global: calculate global stress and strain based on displacements and residual forces of boundary nodes
plot: show graphical representation of mechanical quantities mapped into deformed shape
fem = FE.Model() # create element of class Model
Geometry¶
Define the geometry of the model with method geom. The geometry is subdivided into sections, representing a 1-d or 2-d laminate structure normal to the $x$-direction. The sections are given as a list with the absolute length of each section. The total length of the model is the sum of all section lengths. All lengths are given in units of 1 mm.
Parameters:
sect: list with absolute length of each section in x-direction
LY (optional): height of model in y-direction (default: LY=1)
LZ (optional): thickness of model in z-direction (default: LZ=1)
fem.geom([3, 2, 4]) # define sections in absolute lengths
Materials¶
Material properties are given either in terms of anisotropic elastic constants $C_{11}$, $C_{12}$, $C_{44}$, or – in case of an elastically isotropic material – as Young's modulus $E$ and Poisson's ratio $\nu$. The parameters must be given with their name. The material is created by applying the class Material, which defines a container for the attributes and methods describing material properties and behavior. Elastic properties are defined by calling the method elasticity with the appropriate set of parameters.
Class:
Material
Methods:
elasticity: define linear elastic material either by anisoptropic or by isotropic constants. Note: One complete set of elastic parameters must be provided.
Paramaters:
C11, C12, C44 (optional): anisotropic elastic constants
E, nu (optional): Young's modulus and Poisson's ratio (isotropic elasticity)
one of both parameter sets must be provided
Attributes:
sets internal variables for $C_{11}$, $C_{12}$, $C_{44}$, $E$ and $\nu$
mat1 = Material() # create element of class Material
mat1.elasticity(E=200.e3, nu=0.3) # assign isotropic elastic properties
mat2 = Material()
mat2.elasticity(E=100.e3, nu=0.3)
mat3 = Material()
mat3.elasticity(E=500.e3, nu=0.3)
print('MAT1: C11, C12, C44, E, nu', mat1.C11, mat1.C12, mat1.C44, mat1.E, mat1.nu)
print('MAT2: C11, C12, C44, E, nu', mat2.C11, mat2.C12, mat2.C44, mat2.E, mat2.nu)
print('MAT3: C11, C12, C44, E, nu', mat3.C11, mat3.C12, mat3.C44, mat3.E, mat3.nu)
Assigning materials to model sections¶
Each section has to be assigned a material, which is accomplished by calling the according method of the class Model with the list of materials, which must have the same shape as the list of sections in the geometry definition.
Class: Model
Method:
assign: assign a material to each section<br/assign: assign a material to each section
Paramaters:
mats: list of materials, each defined by class Material, with same shape as in method geom. Each section is assigned to the according material.
Attributes
fem.mat: list with pointer to material for each section
fem.assign([mat1, mat2, mat3]) # assign materials to sections of model
print(fem.mat)
Boundary Conditions¶
For left-hand-side and bottom nodes a displacement is given, typically $u_{\rm left}=u_{\rm bot}=0$.
The type of boundary conditions for right-hand-side and top nodes is specified as either displacement (type='disp') or force (type='force'). Value and type must be specified as parameters in this sequence.
For 1-d models only bcleft and bcright can be called.
Method:
bcleft: set displacement for lhs node(s)
bcbot: set displacement for bottom node(s), has no function for 1-d models
Paramaters:
scalar value: Displacement for boundary node(s)
Method:
bcright: set displacement or force for rhs node(s)
bctop: set displaement or force for top node(s), has no function for 1-d models
Paramaters:
scalar value: Displacement of force for boundary nodes
bctype: string indicating type of boundary conditions: 'disp' for displacements, 'force' for forces
fem.bcleft(0.) # define displacement boundary condition on lhs node (u_x given as parameter)
fem.bcright(0.1*fem.lenx, 'disp') # rhs node is subject to 10% strain (disp = 0.1 * length of model)
Meshing¶
Create a structured mesh with 1-d elements or 2-d quadrilateral elements. Nodes are situated at the corners of the elements and in the middle of the edges for quadratic shape functions. The finite elements of the model are stored in a list with elements of the subclass element with its own attributes and methods, which are used to setup the element shape function.
Method (class Model):
mesh
Parameters:
NX (optional): number of elements in x-direction (default: NX=10)
NY (optional): number of elements in y-direction (default: NY=1), only effective for 2-d models
SF (optional): degree of shape functions: 1=linear (default), 2=quadratic
Subclass: `Element` (parent class: `Model`)
Methods:
calc_Bmat: Calculate B matrix at Gauss point
calc_Kel: calculate element stiffness matrix
Attributes:
nodes: list of nodes of this element
Lelx, Lely: x- and y-dimensions of element
ngp: number of Gauss points
gpx, gpy: np.arrays of x- and y-locations of Gauss points in element
Bmat: list of B-matrices at Gauss points
Vel: volume of element
Kel: element stiffness matrix
Jac: determinant of Jacobian
wght: weight factor for each Gauss point in integration
Sect: Section of model in which element is located
Mat: Pointer to class Material to be applied to this element
CV: Voigt stiffness matrix of element material
eps: average element (total) strain in Voigt notation
sig: average element stress in Voigt notation
fem.mesh(NX=9) # create mesh
print('nodal positions: ', fem.npos.round(decimals=2))
print('nodes belonging to element #0: ', fem.element[0].nodes)
print('stiffness matrix of element #0:')
print(fem.element[0].Kel.round(decimals=2))
Solver¶
The system of linear equation defined by the stiffness matrix is solved, resulting in nodal displacements fulfilling mechanical equilibrium for the given boundary conditions. The nodal forces of internal nodes are zero, forces on boundary nodes are either given as boundary conditions or are residual forces on fixed boundary nodes. If no boundary conditions are specified, the force on such boundary nodes will be zero.
After solution is completed, the following attributes of the class `Model` are defined.
Attributes:
u: list of nodal displacements
f: list of nodal forces
element[i].eps, element[i].sig: Voigt strain and stress tensor of element #i
fem.solve() # solve linear system of equations
print('nodal displacements u = ',fem.u.round(decimals=2))
print('nodal forces f = ',fem.f.round(decimals=2))
print('element strain eps = ',[element.eps[0].round(decimals=2) for element in fem.element])
print('Voigt strain tensor for element #0 = ', fem.element[0].eps.round(decimals=2))
Post-processing¶
Calculate global stresses and strains based on displacements and residual forces on boundary nodes. Provide a graphical output of the specified field, the parameter can be either 'strain' or 'stress'. A rectangle is plotted for each element, with the color given by the specified field and a color map.
Class `Model`
Method:
calc_global: Calculate global stresses and strains based on displacements and residual forces of boundary nodes.
Attributes:
glob: Python dictionary with global strains and stresses, contains the elements:
'ebc1', 'ebc2', 'sbc1', 'sbc2' : global strain and stress calculated from reaction forces and displacements
of boundary nodes (type: float)
'eps', 'epl', 'sig', : global strain, plastic strain, and stress tensors homogenized
from element solutions (type: Voigt tensor)
Method:
plot: Calculate stresses and strains for each element and stores them in lists assigned to model.
Parameters:
fsel: field selector, string with value 'strain1' for $\epsilon_{11}$, 'strain2' for $\epsilon_{22}$, 'stress1' for $\sigma_{11}$, 'stress2' for $\sigma_{22}$.
mag (optional): scaling factor (magnification) for displacements (default: mag=10)
cdepth (optional): depth of colormap (default: cdepth=20)
showmesh (optional): Boolean variable to set/unset plotting of lines for element edges (default: showmesh=True)
shownodes (optional): Boolean variable to set/unset plotting of nodes (default: shownodes=True)
See the online documentation of the pyLabFEA package for a full description of the attributes and methods of the class Model.
fem.calc_global() # calculate global stress and strain
print('global stress: sig_11 =', fem.glob['sig'][0].round(decimals=3), 'MPa')
print('global strain: eps_11 =', (fem.glob['eps'][0]*100).round(decimals=2), '%')
print('\nGraphical output of model with magnification factor 10')
fem.plot('strain1') # create plot
print('\nGraphical output of model with magnification factor 1')
fem.plot('strain1', mag=1)
2-dimensional model 1: iso-strain¶
A 2-d model with three different sections is created and subjected to a uniaxial stress parallel to the sections (y-direction), resulting in an iso-strain condition, because all sections are subject to the same total strain as the entire model. The effective eleastic properties are calculated numerically and compared to the results of the Voigt homogenization rule (iso-strain assumption).
'2-d model, iso-strain'
comp=FE.Model(dim=2, planestress=True)
comp.geom([3,2,4], LY=9) # define sections in absolute lengths
comp.assign([mat1, mat2, mat3]) # assign a material to each section
comp.bcleft(0.) # define boundary conditions
comp.bcright(0., 'force')
comp.bcbot(0.)
comp.bctop(0.1*comp.leny, 'disp') # apply eps_22 strain at top boundary
comp.mesh(NX=18, NY=18) # create structured mesh with 18x18 elements
comp.solve() # solve system of equations for mechanical equilibrium
comp.calc_global() # evaluate global properties
mod_stiff = comp.glob['sig'][1]/comp.glob['eps'][1]
voigt_stiff = 3.*mat1.E/9 + 2*mat2.E/9 + 4*mat3.E/9 # weighted average of Young's moduli wrt volume fractions
print('2-d Model 1: iso-strain, plane stress, eps_22=10%')
print('Global strain: ',comp.glob['eps'][0].round(decimals=3), comp.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', comp.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', comp.element[6*18].eps.round(decimals=3))
print('Element strain Section 3: ', comp.element[10*18].eps.round(decimals=3))
print('Global stress (MPa): ',comp.glob['sig'][0].round(decimals=3), comp.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', comp.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', comp.element[6*18].sig.round(decimals=3))
print('Element stress (MPa) Section 3: ', comp.element[10*18].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', voigt_stiff)
print('Error: ', 1.-mod_stiff/voigt_stiff)
print('----------------------------------------')
comp.plot('strain1', mag=1, shownodes=False)
comp.plot('strain2', mag=1, shownodes=False)
comp.plot('stress1', mag=1, shownodes=False)
comp.plot('stress2', mag=1, shownodes=False)
2-dimensional model 2: iso-stress¶
A 2-d model with three different sections is created and subjected to a uniaxial strain perpendicular to the sections (x-direction) resulting in an iso-stress condition, because all section are subject to the same stress in loading direction as the entire model. The effective eleastic properties are calculated numerically and compared to the results of the Reuss homogenization rule (iso-stress assumption).
'2-d model, iso-stress'
comp2=FE.Model(dim=2, planestress=False)
comp2.geom([3,2,4], LY=9) # define sections in absolute lengths
comp2.assign([mat1, mat2, mat3])
comp2.bcleft(0.)
comp2.bcright(0.1*comp2.lenx, 'disp')
comp2.bcbot(0.)
comp2.bctop(0., 'disp')
comp2.mesh(NX=18, NY=18)
comp2.solve()
comp2.calc_global()
mod_stiff = comp2.glob['sig'][0]/comp2.glob['eps'][0]
reuss_stiff = 1./(3/(9*mat1.C11) + 2/(9*mat2.C11) + 4/(9*mat3.C11)) # Reuss average of stiffness
print('2-d Model 2: isostress, eps_11=10%')
print('Global strain: ',comp2.glob['eps'][0].round(decimals=3), comp2.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', comp2.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', comp2.element[6*18].eps.round(decimals=3))
print('Element strain Section 3: ', comp2.element[10*18].eps.round(decimals=3))
print('Global stress: ',comp2.glob['sig'][0].round(decimals=3), comp2.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', comp2.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', comp2.element[6*18].sig.round(decimals=3))
print('Element stress (MPa) Section 3: ', comp2.element[10*18].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', reuss_stiff)
print('Error: ', 1.-mod_stiff/reuss_stiff)
print('----------------------------------------')
comp2.plot('strain1', mag=1, shownodes=False)
comp2.plot('strain2', mag=1, shownodes=False)
comp2.plot('stress1', mag=1, shownodes=False)
comp2.plot('stress2', mag=1, shownodes=False)
Summary¶
FE models with linear elastic materials have been created and applied to analyze the mechanical response of elastic structures to given boundary conditions. The properties of linear-elastic structures are analysed more closely in the tutorial Composites.
Tutorial 2: Composites¶
In this tutorial, the properties of composites made from different elastic materials are analyzed, and the numerical solution is compared with the expected values from mechanical models.
Micromechanics of Elastic Composites¶
Assess properties of linear-elastic composites by Finite Element Analysis (FEA). Model structures of laminates are created and analysed with the pyLabFEA package, see the online documentation and the tutorial Introduction for detailed information on the functionality of the package. This tutorial uses the matplotlib (https://matplotlib.org/) library for the visualization of results, and NumPy (http://www.numpy.org) for mathematical operations with arrays.
Author: Alexander Hartmaier, ICAMS / Ruhr-Universität Bochum, Germany
March 2020
This work is licensed under a Creative
Commons Attribution-NonCommercial-ShareAlike 4.0 International License
(CC-BY-NC-SA)
The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)
1. Finite Element Model¶
In a first step, the class Model is invoked to generate a container for the Finite Element Model. This container stores all attributes and defines methods to be applied. On this model, the following steps are performed: (i) geometry definition, (ii) material definition, (iii) applying boundary conditions, (iv) meshing, (v) solving of linear equations, and (vi) calculation of global stresses and strains to obtain the effective stiffness of the laminate structure.
Geometry¶
Define a laminate composite geometry for the model with the method geom. The method takes a list with the sizes of the sections of the laminate as input. Later material properties are assigned to each section.
Materials¶
Linear elastic material properties are defined here in terms of Young's modulus $E$ and Poisson ratio $\nu$, by using the class Material. The resulting material definitions are assigned to the sections of the geometry with the assign method of the finite element model.
Boundary Conditions¶
For the 2-d model, different boundary conditions are defined, resulting in longitudinal or transversal loading conditions with respect to the laminate structure. The first case applies a uniaxial stress parallel to the laminate sections, resulting in an iso-strain condition according to the Voigt model. For transversal loading, either uniaxial strain or uniaxial stress perpendicular to the laminate structure is applied. These boundary conditions approximate the iso-stress assumption of the Reuss model to various degrees of fulfillment.
Meshing¶
The laminate model in subdivided into finite elements with nodes at their corners, by calling the method mesh. A mesh with 2-d quadrilateral elements with linear shape functions is created.
Solver¶
By invoking the method solve the system of linear equation defined by the stiffness matrix is solved using the NumPy linalg solver, see http://www.numpy.org. This results in nodal displacements fulfilling mechanical equilibrium for the given boundary conditions, i.e. the total force on any internal node is zero. Forces on boundary nodes are either given as boundary conditions or are residual forces on fixed boundary nodes. If no boundary conditions are specified, the force on that node will be zero.
Postprocessing¶
Calculate and print stresses and strains within elements, revealing the local mechanical condition in each part of the laminate. The graphics are generated with the matplotlib.pyplot libary, see https://matplotlib.org/3.1.0/tutorials/introductory/pyplot.html. Furthermore, the global stresses and strains are calculated, from which the effective stiffness of the laminate is derived.
2. Laminate model under isostrain conditions (Voigt model)¶
In the first step, a laminate model with 5 different sections in generated. The model represents a composite with with a compliant phase (Material: mat1, Young's modulus $E=10$ GPa and volume fraction 75%) and a stiff phase (material: mat2, $E=300$ GPa and volume fraction 25%). This model is subjected to a uniaxial stress longitudinal to the laminate structure to calculate the effective Young's modulus that is compared to the value obtained from the Voigt model.
import pylabfea.model as FE
from pylabfea.material import Material
'laminate model is generated and elastic properties are assigned to each section'
fem_v = FE.Model(dim=2, planestress=True) # call class to generate container for finite element model
fem_v.geom([2, 1, 2, 1, 2], LY=4.) # define sections in absolute lengths
mat1 = Material() # call class to generate container for material
mat1.elasticity(E=10.e3, nu=0.3) # define materials by their elastic propertiess
mat2 = Material() # define second material
mat2.elasticity(E=300.e3, nu=0.3)
fem_v.assign([mat1, mat2, mat1, mat2, mat1]) # assign the proper material to each section
fmat1 = 6./8. # calculate volume fraction of each material
fmat2 = 2./8.
'boundary conditions: uniaxial stress in longitudinal direction'
fem_v.bcleft(0.) # fix left and bottom boundary
fem_v.bcbot(0.)
fem_v.bcright(0., 'force') # free boundary condition on right edge of model
fem_v.bctop(0.1*fem_v.leny, 'disp') # strain applied to top nodes
'solution and evalutation of effective properties'
fem_v.mesh(NX=16, NY=4) # create mesh
fem_v.solve() # solve system of equations
fem_v.calc_global() # calculate global stress and strain
mod_stiff = fem_v.glob['sig'][1]/fem_v.glob['eps'][1] # effective stiffness of numerical model
voigt_stiff = fmat1*mat1.E + fmat2*mat2.E # Voigt stiffness: weighted average of Young's moduli wrt volume fractions
print('2-d Model: isostrain, uniaxial stress, e2=10%')
print('Global strain: ',fem_v.glob['eps'][0].round(decimals=3), fem_v.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', fem_v.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', fem_v.element[16].eps.round(decimals=3))
print('Global stress: ',fem_v.glob['sig'][0].round(decimals=3), fem_v.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', fem_v.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', fem_v.element[16].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', voigt_stiff)
print('Error: ', 1.-mod_stiff/voigt_stiff)
print('----------------------------------------')
fem_v.plot('strain1',mag=1)
fem_v.plot('strain2',mag=1)
fem_v.plot('stress1',mag=1)
fem_v.plot('stress2',mag=1)
It is seen that the isostrain assumption is fulfilled to a very good degree and that the numerical value is identical to that of the analytical Voigt model.
3. Laminate model under isostrain conditions (Reuss model)¶
In the next step, a new model with the identical geometry and materials is generated, but subjected to transversal loading with a uniaxial strain.
'laminate model is generated and elastic properties are assigned to each section'
fem_r = FE.Model(dim=2)
fem_r.geom([2, 1, 2, 1, 2], LY=4.) # define sections in absolute lengths
fem_r.assign([mat1, mat2, mat1, mat2, mat1]) # assign the materials to each section
fmat1 = 6./8. # calculate volume fractions
fmat2 = 2./8.
'mechanical boundary conditions for Reuss model (isostres) with uniaxial strain'
fem_r.bcleft(0.) # fix left and bottom boundary
fem_r.bcbot(0.)
fem_r.bcright(0.1*fem_r.lenx, 'disp') # rhs nodes are subject to 10% strain (disp = 0.1 * length of model)
fem_r.bctop(0., 'disp') # top boundary is fixed
'solution and evaluation of effective properties'
fem_r.mesh(NX=16, NY=4) # create mesh
fem_r.solve() # solve equation for mechanical equilibrium
fem_r.calc_global() # calculate global stress and strain
mod_stiff = fem_r.glob['sig'][0]/fem_r.glob['eps'][0] # calculate effective stiffness of model
reuss_stiff = 1./(fmat1/mat1.C11 + fmat2/mat2.C11) # Reuss average for stiffness
print('2-d Model: isostress, uniaxial strain e1=10%')
print('Global strain: ',fem_r.glob['eps'][0].round(decimals=3), fem_r.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', fem_r.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', fem_r.element[16].eps.round(decimals=3))
print('Global stress: ',fem_r.glob['sig'][0].round(decimals=3), fem_r.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', fem_r.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', fem_r.element[16].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', reuss_stiff)
print('Error: ', 1.-mod_stiff/reuss_stiff)
print('----------------------------------------')
fem_r.plot('strain1',mag=1)
fem_r.plot('strain2',mag=1)
fem_r.plot('stress1',mag=1)
fem_r.plot('stress2',mag=1)
The Reuss model is well fulfilled for uniaxial strain conditions, where the cross-contraction does not change the geometry of the model. Thus the isostress assumption is valid. Under such uniaxial strain conditions, however, the elastic stiffness $C_{11}$ is the relevant quantity and the calculated numerical value is identical to that calculated by the analytical Reuss model. Now we apply a uniaxial stress, to calculate an effective Young's modulus of the model.
'mechanical boundary conditions changed to uniaxial stress'
fem_r.bctop(0., 'force') # top boundary is free (other boundary conditions remain in place)
fem_r.planestress = True # plane stress condition is activated
fem_r.mesh(NX=16, NY=4) # create mesh
'solution and evaluation of effective properties'
fem_r.solve() # solve equation for mechanical equilibrium
fem_r.calc_global() # calculate global stress and strain
mod_stiff = fem_r.glob['sig'][0]/fem_r.glob['eps'][0] # calculate effective Young's modulus of model
reuss_stiff = 1./(fmat1/mat1.E + fmat2/mat2.E) # Reuss average for Young's modulus
print('2-d Model: uniaxial stress, e1=10%')
print('Global strain: ',fem_r.glob['eps'][0].round(decimals=3), fem_r.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', fem_r.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', fem_r.element[16].eps.round(decimals=3))
print('Global stress: ',fem_r.glob['sig'][0].round(decimals=3), fem_r.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', fem_r.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', fem_r.element[16].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', reuss_stiff)
print('Error: ', 1.-mod_stiff/reuss_stiff)
print('----------------------------------------')
fem_r.plot('strain1',mag=1)
fem_r.plot('strain2',mag=1)
fem_r.plot('stress1',mag=1)
fem_r.plot('stress2',mag=1)
Because of the different strains in both phases, their cross contraction is different, which violates the iso-stress assumption. The effective numerical Young's modulus under these conditions deviates from the Reuss value by 5%.
4. Laminate model with different geometrical arrangement of phases¶
In the next step, we compare composites with identical volume fractions, but different geometry, with only two sections, one for each material. Uniaxial stress is applied and the effective Young's modulus in longitudinal and transversal direction is calculated and compared to the Voigt and Reuss model, respectively.
fem = FE.Model(dim=2, planestress=True)
fem.geom([3, 1], LY=4) # create block geometry
fem.assign([mat1, mat2]) # assign materials to sections
f_1 = 3/4 # calculate volume fractions
f_2 = 1/4
'boundary conditions for Voigt model'
fem.bcleft(0.) # fix left and bottom boundary
fem.bcbot(0.)
fem.bcright(0., 'force') # free right edge
fem.bctop(0.1*fem.leny, 'disp') # strain applied to top nodes
'solution and evalutation of effective properties'
fem.mesh(NX=16, NY=4) # create mesh
fem.solve() # calculate solution
fem.calc_global() # calculate global stress and strain
mod_stiff = fem.glob['sig'][1]/fem.glob['eps'][1] # effective Young's modulus of model
voigt_stiff = f_1*mat1.E + f_2*mat2.E # weighted average of Young's moduli wrt volume fractions
print('2-d Model: uniaxial stress, e2=10%')
print('Global strain: ',fem.glob['eps'][0].round(decimals=3), fem.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', fem.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', fem.element[fem.Nel-1].eps.round(decimals=3))
print('Global stress: ',fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', fem.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', fem.element[fem.Nel-1].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', voigt_stiff)
print('Error: ', 1.-mod_stiff/voigt_stiff)
print('----------------------------------------')
fem.plot('stress2',mag=1)
'boundary conditions for Reuss model'
fem.bctop(0., 'force') # free top edge
fem.bcright(0.1*fem.lenx, 'disp') # strain applied to rhs nodes
'solution and evalutation of effective properties'
fem.solve()
fem.calc_global()
mod_stiff = fem.glob['sig'][0]/fem.glob['eps'][0]
reuss_stiff = 1/(f_1/mat1.E + f_2/mat2.E) # weighted average of Young's moduli wrt volume fractions
print('2-d Model: uniaxial stress, e1=10%')
print('Global strain: ',fem.glob['eps'][0].round(decimals=3), fem.glob['eps'][1].round(decimals=3))
print('Element strain Section 1: ', fem.element[0].eps.round(decimals=3))
print('Element strain Section 2: ', fem.element[fem.Nel-1].eps.round(decimals=3))
print('Global stress: ',fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3))
print('Element stress (MPa) Section 1: ', fem.element[0].sig.round(decimals=3))
print('Element stress (MPa) Section 2: ', fem.element[fem.Nel-1].sig.round(decimals=3))
print('Stiffness (MPa): ',mod_stiff)
print('Target (MPa): ', reuss_stiff)
print('Error: ', 1.-mod_stiff/reuss_stiff)
print('----------------------------------------')
fem.plot('strain1',mag=1)
It is seen that the geometrical arragement does not significantly change the results, neither for longitudinal nor for transverse loading.
5. Calculation of effective properties as function of volume fraction¶
The block-like structure is used to calculate the effective elastic properties of the composite as function of the volume fraction.
import numpy as np
import matplotlib.pyplot as plt
'initialize arrays for results'
nmax = 10
f_2 = np.zeros(nmax)
E_long = np.zeros(nmax)
E_trans = np.zeros(nmax)
E_voigt = np.zeros(nmax)
E_reuss = np.zeros(nmax)
for i in range(nmax):
x = (i+1)/(nmax+1)
fem.geom([x, 1.-x])
fem.assign([mat2, mat1])
f_2[i] = x
'boundary conditions for Voigt model'
fem.bcright(0., 'force') # free right edge
fem.bctop(0.1*fem.leny, 'disp') # strain applied to top nodes
'solution and evalutation of effective properties'
fem.mesh(NX=nmax*2, NY=4) # remeshing due to plane stress condition
fem.solve()
fem.calc_global()
E_long[i] = fem.glob['sig'][1]/fem.glob['eps'][1]
E_voigt[i] = (1-x)*mat1.E + x*mat2.E
'boundary conditions for Reuss model'
fem.bcright(0.1*fem.lenx, 'disp') # strain applied to right edge
fem.bctop(0., 'force') # free top nodes
'solution and evalutation of effective properties'
fem.solve()
fem.calc_global()
E_trans[i] = fem.glob['sig'][0]/fem.glob['eps'][0]
E_reuss[i] = 1./((1-x)/mat1.E + x/mat2.E)
'Graphical output'
plt.scatter(f_2,E_long,color='red', label='model_long')
plt.scatter(f_2,E_trans,color='navy', label='model_trans')
plt.plot(f_2,E_voigt,color='darkorange', label='Voigt')
plt.plot(f_2,E_reuss,color='cornflowerblue', label='Reuss')
plt.xlabel('f_2')
plt.ylabel('stiffness (MPa)')
plt.title('Volume fraction vs. stiffness')
plt.legend()
plt.show()
For volume fractions in the range of 10% to 90% the numerical results agree very well with those of the analytical Voigt and Reuss models, even though the iso-stress condition of the Reuss model is not strictly fulfilled.
Tutorial 3: Plasticity¶
Non-linear material behavior in form of plasticity is introduced in this tutorial.
Non-linear FEA for Elastic-Plastic Materials¶
Plastic deformation of materials and structures can be described by non-linear Finite Element Analysis with the pyLabFEA package, see the online documentation and the tutorial Introduction for detailed information on the functionality of the package. This tutorial uses the matplotlib (https://matplotlib.org/) library for the visualization of results.
Author: Alexander Hartmaier, ICAMS / Ruhr-Universität Bochum, Germany
March 2020
This work is licensed under a Creative
Commons Attribution-NonCommercial-ShareAlike 4.0 International License
(CC-BY-NC-SA)
The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)
import pylabfea.model as FE
from pylabfea.material import Material
import matplotlib.pyplot as plt
1. Define plastic material¶
In a first step, the class Material is invoked to define an elastic-plastic material. We start with the easiest case: an isotropic material with ideal plasticity, i.e. no work hardening, is defined. The class Material contains methods to assign elastic and plastic properties to the material, see the online documentation for a complete overview on all attributes and methods of the pyLabFEA package. All stresses are given in units of 1 MPa.
'define elastic-plastic material with isotropic elastic and plastic properties and zero work hardening rate'
E = 200.e3 # standard values for elastic and plastic properties
nu = 0.3
sy = 300.
mat1 = Material(name="isotropic, ideal plastic") # call class to generate object for material
mat1.elasticity(E=E, nu=nu) # define material with isotropic elasticity
mat1.plasticity(sy=sy, khard=0.) # define material with ideal isotropic plasticity
The class material also includes several methods to characterize the material's properties and to generate graphical output. This includes the visualization of the material's yield locus with the method plot_yield_locus, which displays the yield locus as the hyperplane in the principal stress space where the yield function takes the value of zero. The yield function is given as
\begin{equation} f = \sigma_{eq} - \sigma_y\,, \end{equation}
with the initial yield strength of the material $\sigma_y$, and the equivalent stress
which - in this formulation - is based on the pricipal stresses $\sigma_i$ with $i=1, 2, 3$. This equivalent stress follows the definition of von Mises, based on J2, the second invarient of the stress deviator. If the material is exposed to a stress lying within the yield locus, i.e. a stress that produces a negative value of the yield function, the material will respond with a purely elastic deformation. If the stress reaches the yield locus, resulting in a value of zero for the yield function, plastic deformation will set in.
Since we only consider 2D models in this tutorial, the yield locus is drawn in the $\sigma_1$-$\sigma_2$ space, and $\sigma_3=0$.
'Plot yield locus'
mat1.plot_yield_locus(xstart=-1.5, xend=1.5)
It is seen that the typical von Mises-ellipsis is obtained, which is characteristic for isotropic plasticity.
Since the material does not sustain stresses outside the yield locus, plastic flow must occur in a way to keep the stress on the yield locus. Numerically, this is achieved by testing whether an elastic predictor step produces a stress outside the yield locus, i.e., with a positive yield function. If this is the case, a plastic strain increment must be calculated that leads again to an accepted stress state on the yield locus. The return mapping algorithm to calculate such strain increments has been described in many text books on continuum plasticity; the implementation in pyLabFEA follows the book of De Borst, Crisfield, Remmers, and Verhoosel "Nonlinear finite element analysis of solids and structures" (John Wiley & Sons, 2012). According to the Prandtl-Reuss flow rule, the plastic strain increment for a given time step can be calculated as
\begin{equation} \label{eq_epr} \dot\epsilon_{pl} = \dot{\lambda} \frac{\partial f}{\partial \sigma}= \dot{\lambda}{n} \, , \end{equation}where ${n}$ is the normal vector to the yield locus, defined by the gradient of the yield function $\partial f/\partial \sigma_i$ ($i=1,2,3)$, and $\dot{\lambda}>0$ is the so-called plastic strain multiplier that can be evaluated as
\begin{equation} \dot{\lambda} = \frac{{n} \cdot {C}{\dot{\epsilon}}}{{n} \cdot {C} {n}} \, , \end{equation}where $\dot{\epsilon}$ is the total strain increment of the FEA predictor step that leads to a stress state outside the yield locus and which is consequently decomposed into the plastic strain increment and the elastic strain increment - or stress increment - given by
\begin{equation} \dot{\sigma} = {C}_t \dot{\epsilon} \end{equation}with the tangent stiffness tensor
\begin{equation} {C}_t = {C} - \frac{{C} {n} \otimes {C} {n}} {{n} \cdot {C} {n}} \end{equation}where '$\otimes$' denotes the tensorial product in the form $a_i \otimes b_j = a_i b_j$.
The method calc_properties can be invoked to evaluate the flow behavior of the defined material under different load cases: (i) uniaxial stress in horizontal direction, (ii) uniaxial stress in vertical direction, (iii) equibiaxial strain under plane stress conditions, and (iv) pure shear under plane stress conditions. This is achieved by internally calling a finite element simulation to apply these load cases to the material. To plot the obtained stress-strain curves, the method plot_stress_strain is called, which also produces a text output.
'Calculate and plot stress strain curves of materials'
mat1.calc_properties(eps=0.05, sigeps=True) # calculate the stress-strain curves up to a total strain of 5%
mat1.plot_stress_strain()
As expected for an isotropic material, all load cases result in the same stress-strain curve, when plotted as equivalent stress vs. equivalent strain.
2. Work hardening¶
In the next step, two more materials will be defined, with identical elastic properties and identical yield strength as the first material, but with different work hardening behavior, as seen in the stress-strain curves.
'Define two more materials with different work-hardening rates'
mat2 = Material(name="low work-hardening") # define second material
mat2.elasticity(E=200.e3, nu=0.3) # identic elastic properties as mat1
mat2.plasticity(sy=300., khard=10.e3) # ssame yield strength as mat1, low w.h. coefficient
mat3 = Material(name="strong work-hardening") # define third material
mat3.elasticity(E=200.e3, nu=0.3) # identic elastic properties as mat1
mat3.plasticity(sy=300., khard=30.e3) # same yield strength as mat1 and mat2, high w.h. coefficient
mat2.calc_properties(verb=False, eps=0.05, sigeps=True)
mat3.calc_properties(verb=False, eps=0.05, sigeps=True)
plt.plot(mat1.prop['stx']['eeq']*100., mat1.prop['stx']['seq'], '-b')
plt.plot(mat2.prop['stx']['eeq']*100., mat2.prop['stx']['seq'], '-m')
plt.plot(mat3.prop['stx']['eeq']*100., mat3.prop['stx']['seq'], '-r')
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon$ (%)')
plt.ylabel(r'$\sigma$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name], loc='lower right')
plt.show()
3. Application in Finite Element Analysis¶
The materials defined are used in a finite element model, in which 3 different sections are assigned to the 3 different materials. This model is subjected to a unixial tensile stress in vertical direction, and the mechanical response of the model is analyzed in terms of total strain, plastic strain, and stress. All lengths are given in units of 1 mm.
'laminate model is generated and elastic-plastic properties are assigned to each section'
fem = FE.Model(dim=2, planestress=True) # call class to generate object for finite element model
fem.geom([2, 2, 2], LY=6.) # define sections in absolute lengths
fem.assign([mat1, mat2, mat3]) # assign the proper material to each section
fmat1 = 1./3. # identical volume fraction for each material
fmat2 = 1./3.
fmat3 = 1./3.
'boundary conditions: uniaxial stress in longitudinal direction'
etot = 0.002
fem.bcleft(0.) # fix left and bottom boundary
fem.bcbot(0.)
fem.bcright(0., 'force') # free boundary condition on right edge of model
fem.bctop(etot*fem.leny, 'disp') # strain applied to top nodes
'solution and evalutation of effective properties'
fem.mesh(NX=6, NY=4) # create mesh
fem.solve() # solve system of equations
fem.calc_global() # calculate global stress and strain
print('2-d Model: uniaxial stress, isostrain condition')
print('-----------------------------------------------')
print('Global total strain: (eps_xx, eps_yy) = (%6.3f, %6.3f)' % (fem.glob['eps'][0], fem.glob['eps'][1]))
print('Local total strain (elements in different sections, Voigt tensor):')
print(fem.element[0].eps.round(decimals=3),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[8].eps.round(decimals=3),' in Section 2, Material:',fem.element[8].Mat.name)
print(fem.element[16].eps.round(decimals=3),' in Section 3, Material:',fem.element[16].Mat.name)
'create graphical output of total strains in FE model'
print('Total strain in vertical direction: isostrain condition fulfilled')
fem.plot('strain2')
print('Total strain in horizontal direction:')
print('slightly different cross-contraction due to different plastic properties in model sections')
fem.plot('strain1')
print('Equivalent total strain')
fem.plot('etot')
Since isostrain conditions are applied, the global total strain in vertical direction is identical in all elements, irrespective of their section and the assigned material. However, due to the differences in work hardening, slight changes in the horizontal strain show up. Next, we take a look at the plastic strains:
print('Global plastic strain: (epl_xx, epl_yy) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0], fem.glob['epl'][1]))
print('Local plastic strain (elements in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 2, Material:',fem.element[8].Mat.name)
print(fem.element[16].epl.round(decimals=6),' in Section 3, Material:',fem.element[16].Mat.name)
'create graphical output of plastic strains in FE model'
print('Plastic strain in vertical direction: differences due to work hardening')
fem.plot('plastic2',vmin=0.047,vmax=0.051)
print('Plastic strain in horizontal direction: different plastic cross-contraction')
fem.plot('plastic1')
print('Equivalent plastic strain')
fem.plot('peeq',vmin=0.047,vmax=0.051)
print('Global stress: (sig_xx, sig_yy) = (%6.3f, %6.3f) MPa'
% (fem.glob['sig'][0], fem.glob['sig'][1]))
print('Local stress (elements in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=3),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[8].sig.round(decimals=3),' in Section 2, Material:',fem.element[8].Mat.name)
print(fem.element[16].sig.round(decimals=3),' in Section 3, Material:',fem.element[16].Mat.name)
'create graphical output of stress in FE model'
print('Stress in vertical direction: different flow stresses at constant strain')
fem.plot('stress2',vmin=300,vmax=305)
print('Stress in horizontal direction: uniaxial stress condition fulfilled')
fem.plot('stress1')
print('Equivalent stress: different flow stresses at constant strain')
fem.plot('seq',vmin=300,vmax=305)
4. Anisotropic materials¶
The pyLabFEA package also allows the definition of plastically anisotropic materials. When defining the plastic properties, the parameter hill$=[H_1, H_2, H_3]$ can be used to define orthotropic plastic properties in a Hill-like manner. The parameters are taken into account in a generalized equivalent stress
In the first step, the yield locus of an anisotropic material is plotted together with the isotropic reference material mat1. Furthermore, the stresses obtained during the plastic deformation in the internal calculation of the stress-strain curves under the different loading conditions are also shown.
'Material with anistropic plasticity as described by Hill-parameters'
tol = FE.ptol # save tolerance original parameter
FE.ptol = 0.4 # increase tolerance value to speed up calculation for strongly anisotropic material
mat_an = Material(name='anisotropic Hill')
mat_an.elasticity(E=E, nu=nu)
mat_an.plasticity(sy=sy, hill=[0.7,1.,1.4])
mat_an.calc_properties(eps=0.005, sigeps=True)
FE.ptol = tol # reset tolerance to its original value
'plot yield locus'
ax = mat_an.plot_yield_locus(xstart=-1.5, xend=1.5, ref_mat=mat1)
'plot evolution of stresses during loading for anisoptropy ...'
stx = mat_an.sigeps['stx']['sig'][:,0:2]/mat_an.sy
sty = mat_an.sigeps['sty']['sig'][:,0:2]/mat_an.sy
et2 = mat_an.sigeps['et2']['sig'][:,0:2]/mat_an.sy
ect = mat_an.sigeps['ect']['sig'][:,0:2]/mat_an.sy
ax.scatter(stx[1:,0],stx[1:,1],s=60, c='#00f0ff', edgecolors='k')
ax.scatter(sty[1:,0],sty[1:,1],s=60, c='#00f0ff', edgecolors='k')
ax.scatter(et2[1:,0],et2[1:,1],s=60, c='#0000ff', edgecolors='k')
ax.scatter(ect[1:,0],ect[1:,1],s=60, c='#0060ff', edgecolors='k')
'... and for the isotropic material mat1'
stx = mat1.sigeps['stx']['sig'][:,0:2]/mat1.sy
sty = mat1.sigeps['sty']['sig'][:,0:2]/mat1.sy
et2 = mat1.sigeps['et2']['sig'][:,0:2]/mat1.sy
ect = mat1.sigeps['ect']['sig'][:,0:2]/mat1.sy
ax.scatter(stx[1:,0],stx[1:,1],s=60, c='#00f0ff', edgecolors='k')
ax.scatter(sty[1:,0],sty[1:,1],s=60, c='#00f0ff', edgecolors='k')
ax.scatter(et2[1:,0],et2[1:,1],s=60, c='#ffff00', edgecolors='k')
ax.scatter(ect[1:,0],ect[1:,1],s=60, c='#ffff00', edgecolors='k')
It is seen how the shape of the yield locus gets deformed by the anisotropy paramaters. Furthermore, the graph shows that all obtained stress values during the calculation of the stress-strain curves lie on the yield locus. For both materials shown here, ideal plasticity is assumed such that the yield locus is not changed during plastic deformation. For the anisotropic material, the flow stresses obtained under equibiaxial strain and pure shear strain move on the yield locus, which is a consequence of the material anisotropy that changes the stress state during strain controled loading conditions. For loading cases leading to uniaxial stresses, the flow stress remains on on the same spot of the yield locus by definition.
Finally, the obtained stress strain courves are plotted for the two definitions of the equivalent stress: (i) isotropic J2 equivalent stress, (ii) Hill-like anisotropic equivalent stress.
mat_an.plot_stress_strain(Hill=True)
It can be seen that for anisotropy the two different ways to define the equivalent stress, result in different values describing the same material under the same loading conditions. While the J2 equivalent stress reveales the anisotropy of the material in different yield stresses, depending on the loading direction, the equivalent stress after Hill maps the yield stresses on a constant strength value. Note, that the nominal yield strength of the material is the constant in both cases.
5. Summary¶
This tutorial demonstrates the non-linear FEA with the pyLabFEA package. Yield loci and stress-strain curves for materials with isotropic plasticity and Hill-like anisotropic behavior are analysed. FE models with materials exhibiting different work hardening rates are generated and applied.
6. Exercises¶
- Create materials with different initial strengths and different work hardening rates and compare the results of the FE model.
- Create materials with identical plastic properties but different elasticity constants and compare the results of the FEA for loading in different directions.
- Combine different materials with identical yield strengths and work hardening rates, but differnt Hill parameters in one FE model and discuss the results.
Tutorial 4: Homogenization¶
Laminate structures with different elastic-plastic materials are analyzed with respect to the mechanical behavior.
Homogenization of non-linear materials¶
The mechanical properties of structures composed of materials with different elasto-plastic properties can be estimated based on the models of Taylor or Sachs. In this tutorial, structures composed of different materials are created and their mechanical properties are investigated and compared to the Taylor and Sachs solution. Non-linear Finite Element Analysis is performed with the pyLabFEA package, see the online documentation and the tutorial Introduction for detailed information on the functionality of the package. For linear-elastic materials, analytical solutions for the homogenzition of their properties exist, as described in the tutorial Composites. This tutorial uses the matplotlib (https://matplotlib.org/) library for the visualization of results, and NumPy (http://www.numpy.org) for mathematical operations with arrays.
Author: Alexander Hartmaier, ICAMS / Ruhr-Universität Bochum, Germany
March 2020
This work is licensed under a Creative
Commons Attribution-NonCommercial-ShareAlike 4.0 International License
(CC-BY-NC-SA)
The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)
import pylabfea.model as FE
from pylabfea.material import Material
import matplotlib.pyplot as plt
import numpy as np
1. Create elastic-plastic materials¶
In the first step, the class Material in invoked to create objects for three materials, with identical elastic properties, but different yield strengths. The unit of stress is 1 MPa. All features of the class Material are described in the documentation.
mat1 = Material(name="Soft") # invoke class to generate object for material 1
mat1.elasticity(E=200.e3, nu=0.3) # define material with isotropic elasticity
mat1.plasticity(sy=180., khard=50.e3) # define material with isotropic plasticity and linear work hardening
mat2 = Material(name="Intermediate") # define second material
mat2.elasticity(E=200.e3, nu=0.3) # identic elastic properties as mat1
mat2.plasticity(sy=190., khard=50.e3) # higher yield strength than mat1, but same w.h. coefficient
mat3 = Material(name="Strong") # define third material
mat3.elasticity(E=200.e3, nu=0.3) # identic elastic properties as mat1
mat3.plasticity(sy=230., khard=50.e3) # higher yield strength than mat1 and mat2, but same w.h. coefficient
'Calculate and plot stress strain curves of materials'
mat1.calc_properties(eps=0.01) # invoke method to calculate properties of material
mat2.calc_properties(eps=0.01) # and to store stress-strain data up to total strain eps=1%
mat3.calc_properties(eps=0.01)
plt.plot(mat1.prop['sty']['eeq']*100., mat1.prop['stx']['seq'], '-b') # plot equiv. strain (eeq) vs. equiv. stress (seq)
plt.plot(mat2.prop['sty']['eeq']*100., mat2.prop['stx']['seq'], '-m') # for uniaxial tensile stress in vertical direction (sty)
plt.plot(mat3.prop['sty']['eeq']*100., mat3.prop['stx']['seq'], '-r')
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{eq}$ (%)')
plt.ylabel(r'$\sigma_{eq}$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name], loc='lower right')
plt.show()
2. Laminate structure under iso-strain condition (Taylor model)¶
The class Model is invoked to generate an object for FEA. The methods of the class are used to generate the geomtry, to assign materials to sections, to establish boundary conditions, and to create the mesh. After these pre-processing steps, the solver is called to determine the solution for the deformed shape of the laminate structure in mechanical equilibrium under the applied boundary conditions. In the post-processing step, the results are evaluated and visualized. See the tutorial Introduction for the basic steps; all features of the class Model are described in the documentation.
The FE model of the laminate structure is subjected to a uniaxial stress along the lamellae of the structure, leading to iso-strain conditions, in which each section of the model, i.e. each lamella, undergoes the same total strain as the entire structure. The unit of length is 1 mm. In the first step, the model is loaded up to a stress slightly higher than the yield strength of mat1, and the strains are analyzed.
'laminate model is generated and elastic-plastic properties are assigned to each section'
fem = FE.Model(dim=2, planestress=True) # generate an object for a 2d finite element model with plane stress conditions
fem.geom([2, 2, 2], LY=4.) # define 3 horizontal sections with 2 micron width per section, vertical dimension is 4 micron
fem.assign([mat1, mat2, mat3]) # assign the proper material to each section
fmat1 = 1./3. # identical volume fraction for each material
fmat2 = 1./3.
fmat3 = 1./3.
'boundary conditions: uniaxial stress along the direction of the lamellae'
fem.bcleft(0.) # fix left and bottom boundary in normal directions
fem.bcbot(0.) # nodes on lhs boundary can move up and down, bottom nodes can move left and right
fem.bcright(0., 'force') # free boundary condition on right edge of plane stress model -> uniaxial stress condition
ubc = 1.05*fem.leny*mat1.sy/mat1.E # displacement on top boundary corresponding to yield strength of mat1 (+5%)
fem.bctop(ubc, 'disp') # displacement applied to top nodes
'generate mesh'
fem.mesh(NX=6, NY=2) # create structured mesh
'solution of equations for mechanical equilibrium'
fem.solve() # solve system of equations
'post-processing'
fem.calc_global() # calculate global stress and strain
print('2-d Model: iso-strain, uniaxial stress')
print('--------------------------------------')
print('Global total strain: (eps_11, eps_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
'create graphical output of model'
fem.plot('strain2') # plot eps_22
fem.plot('peeq',vmin=0.,vmax=0.0015) # plot eqiv. plastic strain
It is seen that the isostrain assumption is fulfilled to a very good degree, although the first section of the model, to which the soft material mat1 is assigned starts to yield plastically, while the other sections remain linear-elastic.
The analysis of the stresses reveals that the stress in section 1 remains at the level of the flow stress in mat1 (=yield strength + work hardening), while the stress in the two other sections is below the yield strength of those materials.
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('seq')
In the next step the load is increased to exceed the yield strength of mat2, and the strains are analysed again. Note that only the boundary conditions need to be changed on the existing model.
'increase load up to yield strength of mat2'
ubc = 1.05*fem.leny*mat2.sy/mat1.E # displacement on top boundary corresponding to yield strength of mat2 (+5%)
fem.bctop(ubc, 'disp') # displacement applied to top nodes
fem.solve() # solve system of equations
fem.calc_global() # calculate global stress and strain
print('2-d Model: isostrain, uniaxial stress')
print('-------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
'create graphical output of model'
fem.plot('strain2')
fem.plot('peeq',vmin=0.,vmax=0.003)
Again, the iso-strain assumption is fulfilled very well, because all sections individually are subject to the same vertical strain as the entire structure. The plastic strain in section 1 increases, and also in section 2 plastic yielding sets in.
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('seq')
The stress in sections 1 and 2 is bounded by the flow stress of the assigned materials, the stress in section 3 is still below the yield strength of mat3.
Finally, the strain applied on the model is increased to $\epsilon_{22} = 1$%.
'increase load up to eps_22=1%'
ubc = 0.01*fem.leny # displacement on top boundary corresponding to 1% total strain
fem.bctop(ubc, 'disp') # displacement applied to top nodes
fem.solve() # solve system of equations
fem.calc_global() # calculate global stress and strain
print('2-d Model: isostrain, uniaxial stress')
print('-------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
'create graphical output of model'
fem.plot('strain2')
fem.plot('peeq',vmin=0.,vmax=0.008)
The deformation of the laminate structure is now becoming visible, and the plastic strains in the different sections are rather similar, because most of the deformation is now of plastic nature, whereas the elastic contributions is quite small.
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('seq')
The yield strength is now exceed in all sections, due to the linear work hardening with identical rate for each material, the difference in the yield strength is still reflected in the equivalent stress in each section, which corresponds to the flow stress of the assigned materials.
In the next step, the equiv. stress vs. equiv. strain diagram for the laminate structure is plotted together with those of the three materials.
'plot stress-strain curves of materials and laminate structure'
plt.plot(mat1.prop['sty']['eeq']*100., mat1.prop['stx']['seq'], '-b')
plt.plot(mat2.prop['sty']['eeq']*100., mat2.prop['stx']['seq'], '-m')
plt.plot(mat3.prop['sty']['eeq']*100., mat3.prop['stx']['seq'], '-r')
plt.plot(FE.eps_eq(fem.egl)*100., FE.seq_J2(fem.sgl), '-k') # calculate and plot equivalent values of global strain and stress
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{eq}$ (%)')
plt.ylabel(r'$\sigma_{eq}$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name, 'laminate structure'], loc='lower right')
plt.show()
It is seen that the global stress for the laminate structure corresponds to the average of the stresses in the three sections, given by the flow stress of the assigned materials. According to the Taylor model, the stress average is calculated at a constant total strain $\epsilon_\mathrm{tot}$, which is equal to the total strain $\epsilon^{(i)}$ in each section $i$ (iso-strain assumption), thus
\begin{equation} \epsilon^{(i)} = \epsilon_\mathrm{tot} \hspace{2em} (i=1,2,3). \end{equation}The global stress $\sigma_0$ is obtained by averaging the stresses $\sigma^{(i)}$ of the section according to their volume fractions $f_i$ , which yields
\begin{equation} \sigma_0 = \sum\limits_{i=1}^3 f_i \, \sigma^{(i)}(\epsilon_\mathrm{tot}). \end{equation}In the next cell, the stresses according to the Taylor model are calculated and plotted together with the FE solution. To accomplish this, the FE model is restarted an the element solution in each section is recorded for small load increments until the total strain is reached.
'initialize arrays for element stresses and strains'
sig1=[]
sig2=[]
sig3=[]
etot=[]
'reset model to start with zero load; geometry, sections, boundary conditions and mesh are kept'
fem.u=None # reset model
'calculate FE solution for small load increments'
for eps in [8.e-4, 1.e-3, 1.3e-3, 2.e-3, 4.e-3, 6.e-3, 8.e-3, 1.e-2]:
ubc = eps*fem.leny # calculate load increment
fem.bctop(ubc, 'disp') # apply displacement on top nodes
fem.solve() # solve system of equations
sig1.append(fem.element[0].sig[1]) # store element solution for each section
sig2.append(fem.element[4].sig[1])
sig3.append(fem.element[8].sig[1])
etot.append(fem.element[0].eps[1]*100)
'calculate weigthed average of stresses in different sections at constant total strains'
s_taylor = fmat1*np.array(sig1) + fmat2*np.array(sig2) + fmat3*np.array(sig3)
'plot stresses and strains'
plt.plot(etot, sig1, '-b')
plt.plot(etot, sig2, '-m')
plt.plot(etot, sig3, '-r')
plt.plot(fem.egl[:,1]*100., fem.sgl[:,1], '-k') # plot vertical components of global strain and stress
plt.plot(etot, s_taylor, 'oc')
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{22}$ (%)')
plt.ylabel(r'$\sigma_{22}$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name, 'FE solution', 'Taylor model'], loc='lower right')
plt.show()
It is seen that the results of the Taylor model and the FE solution are in excellent agreement. Furthermore, the element solution for the uniaxial stress state corresponds to the stress-strain curves obtained for the individual materials.
3. Laminate model under iso-stress conditions (Sachs model)¶
In the next step, a new model with the identical geometry and materials is generated, but subjected to a uniaxial strain in the direction transversal to lamellae. This load case results in an iso-stress condition that can be described with the Sachs model.
'mechanical boundary conditions for Sachs model (iso-stress) with uniaxial strain'
fem = FE.Model(dim=2, planestress=False) # generate an object for a 2d finite element model with plane strain conditions
fem.geom([2, 2, 2], LY=4.) # define 3 horizontal sections with 2 micron width per section, vertical dimension is 4 micron
fem.assign([mat1, mat2, mat3]) # assign the proper material to each section
fmat1 = 1./3. # identical volume fraction for each material
fmat2 = 1./3.
fmat3 = 1./3.
'boundary conditions: uniaxial stress in transversal direction'
fem.bcleft(0.) # fix left and bottom boundary in normal directions
fem.bcbot(0.) # nodes on lhs boundary can move up and down, bottom nodes can move left and right
ubc = 0.0015*fem.lenx # apply uniaxial strain eps_11 = 0.15%
fem.bcright(ubc, 'disp') # apply displacement on rhs nodes
fem.bctop(0., 'disp') # top boundary is fixed for uniaxial strain
'generate mesh'
fem.mesh(NX=6, NY=2) # create structured mesh
'solution of equations for mechanical equilibrium'
fem.solve() # solve system of equations
'post-processing'
fem.calc_global() # calculate global stress and strain
print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('strain1')
fem.plot('strain2')
fem.plot('peeq',vmin=0., vmax=0.008)
It is seen that the uniaxial strain condition is applied and that sections 1 and 2 already underwent plastic deformation.
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('stress1')
fem.plot('seq')
Also the iso-stress condition in loading direction is fullfilled. Transversal stresses occur due to the restricted cross-contraction, they are different in each section, due to their different plastic deformation, which is reflected in the different equivalent stresses in each section. Note also, that the stress state has a significant triaxiality, because the stress component in x-direction is larger than the equivalent stress, which is independent of hydrostatic stress components.
Now, the strain is increased to $\epsilon_\mathrm{tot} = 1%$, and the results are analyzed.
'increase load'
ubc = 0.01*fem.lenx
fem.bcright(ubc, 'disp') # displacement on rhs nodes
fem.solve() # solve equations for mechanical equilibrium
fem.calc_global() # calculate global stress and strain
print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('strain1')
fem.plot('peeq',vmin=0., vmax=0.45)
Similar to the iso-strain conditions, the model is now fully plastified. Under uniaxial strain, however, the plastic deformation is smaller than for uniaxial stress, due to the pronounced stress triaxiality that is caused by the restricted cross-contraction.
Looking at the stresses reveals these strong components transversal to the loading direction.
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('stress1')
fem.plot('stress2')
fem.plot('seq')
The iso-stress condition in loading direction is fullfilled, i.e. the horizontal stress component is constant throughout the model. However, it is seen that due to the severe stress triaxiality, there are different equivalent stresses in the different sections, following the strengths of the assigned materials.
The homogenization under iso-stress conditions is covered by the Sachs model, where the stress in each section $\sigma^{(i)}$ equals the global stress $\sigma_0$, thus
\begin{equation} \sigma^{(i)} = \sigma_0 \hspace{2em} (i=1,2,3). \end{equation}The global total strain is given as the weighted average of the total strain in each section $\epsilon_\mathrm{tot}^{(i)}$ at a constant stress $\sigma_0$ with the volume fractions $f_i$ as weights, which yields
\begin{equation} \epsilon_\mathrm{tot} = \sum\limits_{i=1}^3 f_i \, \epsilon_\mathrm{tot}^{(i)}(\sigma_0). \end{equation}In the next cell, the strains according to the Sachs model are calculated from the element solutions in each section and plotted together with the global solution.
sig0=[]
seq1=[]
seq2=[]
seq3=[]
et1=[]
et2=[]
et3=[]
eq1=[]
eq2=[]
eq3=[]
fem.u=None # reset model
for eps in [1.e-3, 1.2e-3, 1.4e-3, 1.6e-3, 2.e-3, 3.e-3, 6.5e-3, 1.e-2]:
ubc = eps*fem.lenx # uniaxial strain increments
fem.bcright(ubc, 'disp') # apply displacement on rhs nodes
fem.solve() # solve system of equations
sig0.append(fem.element[0].sig[0])
et1.append(fem.element[0].eps[0]*100)
et2.append(fem.element[4].eps[0]*100)
et3.append(fem.element[8].eps[0]*100)
seq1.append(FE.seq_J2(fem.element[0].sig))
seq2.append(FE.seq_J2(fem.element[4].sig))
seq3.append(FE.seq_J2(fem.element[8].sig))
eq1.append(FE.eps_eq(fem.element[0].eps)*100)
eq2.append(FE.eps_eq(fem.element[4].eps)*100)
eq3.append(FE.eps_eq(fem.element[8].eps)*100)
'weigthed average of strains in different sections at same stress level'
e_sachs = fmat1*np.array(et1) + fmat2*np.array(et2) + fmat3*np.array(et3)
'plot stresses and strains'
plt.plot(et1, sig0, '-b')
plt.plot(et2, sig0, '-m')
plt.plot(et3, sig0, '-r')
plt.plot(fem.egl[:,0]*100, fem.sgl[:,0], '-k') # plot horizontal components of global strain and stress
plt.plot(e_sachs, sig0, 'oc')
plt.xlim((0,0.3))
plt.ylim((0,700))
plt.show()
The numerical solution and the Sachs model are in excellent agreement. Due to the large stress triaxility, the values of the flow stress are quite high, such that the differences between the different sections, caused by the comparatively small differences in the yield strength of each material, appear to be small in comparison. Thus, the evolution of equivalent stresses and equivalent strains in the elemnts of each section is also evaluated.
'weigthed average of equiv. strains and stresses in different sections'
eeq_sachs = fmat1*np.array(eq1) + fmat2*np.array(eq2) + fmat3*np.array(eq3)
seq_sachs = fmat1*np.array(seq1) + fmat2*np.array(seq2) + fmat3*np.array(seq3)
'plot stresses and strains'
plt.plot(eq1, seq1, '-b')
plt.plot(eq2, seq2, '-m')
plt.plot(eq3, seq3, '-r')
plt.plot(FE.eps_eq(fem.egl)*100., FE.seq_J2(fem.sgl), '-k') # plot horizontal components of global strain and stress
plt.plot(eeq_sachs, seq_sachs, 'oc')
plt.title('Stress-strain curves')
plt.xlabel(r'$\epsilon_{11}$ (%)')
plt.ylabel(r'$\sigma_{11}$ (MPa)')
plt.legend([mat1.name, mat2.name, mat3.name, 'FE solution', 'Taylor model'], loc='lower right')
plt.show()
The analysis of eqivalent stresses and strains makes the influence of the material properties more visible. In the next plot, Taylor and Sachs solutions are compared.
plt.plot(etot, s_taylor, '-b')
plt.plot(eeq_sachs, seq_sachs, '-r')
plt.title('Taylor vs. Sachs model')
plt.xlabel(r'$\epsilon_\mathrm{tot}$ (%)')
plt.ylabel(r'$\sigma_0$ (MPa)')
plt.legend(['Taylor model', 'Sachs model'], loc='lower right')
plt.ylim((0,350))
plt.show()
The Taylor model (iso-strain) results in a stronger material behavior than the Sachs model (iso-stress), in which the softer regions dominate the behavior of the entire structure.
Iso-stress approximation by uniaxial stress¶
In the next step, we investigate how applying a uni-axial stress transversal to the laminate structure changes the results.
'laminate model is generated and elastic-plastic properties are assigned to each section'
fem = FE.Model(dim=2, planestress=True) # generate an object for a 2d finite element model with plane stress conditions
fem.geom([2, 2, 2], LY=4.) # define 3 horizontal sections with 2 micron width per section, vertical dimension is 4 micron
fem.assign([mat1, mat2, mat3]) # assign the proper material to each section
fmat1 = 1./3. # identical volume fraction for each material
fmat2 = 1./3.
fmat3 = 1./3.
'boundary conditions: uniaxial stress in transversal direction'
fem.bcleft(0.) # fix left and bottom boundary in normal directions
fem.bcbot(0.) # nodes on lhs boundary can move up and down, bottom nodes can move left and right
fem.bctop(0., 'force') # free boundary condition on top edge of plane stress model -> uniaxial stress condition
ubc = 1.05*fem.lenx*mat1.sy/mat1.E # displacement on rhs boundary corresponding to yield strength of mat1 (+5%)
fem.bcright(ubc, 'disp') # displacement applied to rhs nodes
'generate mesh'
fem.mesh(NX=6, NY=4) # create structured mesh with higer number of elements in verical direction
'solution of equations for mechanical equilibrium'
fem.solve() # solve system of equations
'post-processing'
fem.calc_global() # calculate global stress and strain
print('2-d Model: transversal uniaxial stress')
print('--------------------------------------')
print('Global total strain: (eps_11, eps_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
'create graphical output of model'
fem.plot('strain1') # plot eps_11
fem.plot('strain2') # plot eps_22
fem.plot('peeq',vmin=0.,vmax=0.0015) # plot eqiv. plastic strain
Plastic yielding starts again in section 1, in which also the total strain localises. This causes also a slight difference in the cross contraction of section 1 compared to the other sections.
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('stress1')
fem.plot('stress2')
The different cross-contractions of the sections cause significant stress components transversal to the loading direction. hence, altough the applied boundary condition correspond to a uniaxial stress, within the structure, the stress state deviates from the applied one.
'increase load above the yield strength of mat2'
ubc = 1.3*fem.lenx*mat2.sy/mat2.E
fem.bcright(ubc, 'disp') # displacement on rhs nodes
fem.solve() # solve equations for mechanical equilibrium
fem.calc_global() # calculate global stress and strain
print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('peeq',vmin=0., vmax=0.015)
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('stress1')
Upon further loading, the work hardening in section 1 causes the flow stress in this region to increase. Finally, the stress in the laminate structure exceeds the yield strength of the material in section 2. The iso-stress condition is fulfilled to some approximation in the structure, although the different sections reveal pronounced differences in the cross-contraction.
'increase load to eps_11=1%'
ubc = 0.01*fem.lenx
fem.bcright(ubc, 'disp') # displacement on rhs nodes
fem.solve() # solve equations for mechanical equilibrium
fem.calc_global() # calculate global stress and strain
print('2-d Model: iso-stress, uniaxial strain')
print('--------------------------------------')
print('Global total strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['eps'][0].round(decimals=4), fem.glob['eps'][1].round(decimals=4)))
print('Local total strain (element solution in different sections, Voigt tensor)')
print(fem.element[0].eps.round(decimals=4),'in Section 1, Material: ', fem.element[0].Mat.name)
print(fem.element[4].eps.round(decimals=5),'in Section 2, Material: ', fem.element[4].Mat.name)
print(fem.element[8].eps.round(decimals=5),'in Section 3, Material: ', fem.element[8].Mat.name)
print('\nGlobal plastic strain: (epl_11, epl_22) = (%9.6f, %9.6f)'
% (fem.glob['epl'][0].round(decimals=3), fem.glob['epl'][1].round(decimals=3)))
print('Local plastic strain (element solution in different sections, Voigt tensor):')
print(fem.element[0].epl.round(decimals=6),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].epl.round(decimals=6),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].epl.round(decimals=6),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('strain1')
fem.plot('peeq',vmin=0., vmax=0.85)
At the maximum strain, the laminate structure is again fully plastified. The total strain in loading direction reflects the different states of work hardening in the materials assigned to the different sections.
print('\nGlobal stress (MPa): (sig_11, sig_22) = (%5.2f, %5.2f)'
% (fem.glob['sig'][0].round(decimals=3), fem.glob['sig'][1].round(decimals=3)))
print('Local stress (element solution in different sections, Voigt tensor) (MPa):')
print(fem.element[0].sig.round(decimals=2),' in Section 1, Material:',fem.element[0].Mat.name)
print(fem.element[4].sig.round(decimals=2),' in Section 2, Material:',fem.element[4].Mat.name)
print(fem.element[8].sig.round(decimals=2),' in Section 3, Material:',fem.element[8].Mat.name)
fem.plot('stress1')
The iso-stress assumption is still fulfilled within an error margin of 5%, despite the different cross contractions
4. Summary¶
This tutorial demonstrates the homogenization of structures composed of materials with different elastic-plastic properties by non-linear FEA. Models fhe different load cases leading to iso-strain and iso-stress conditions in the structures are set up and the results are analyzed. It is seen that the numerical solution agrees very well with the analyitcal homogenization rules according to the Taylor and Sachs models.
5. Exercises¶
- Create materials with different hardening rates, e.g. 100 GPa for
mat1, 50 GPa format2and 10 GPa format3and observe how the homogenized stress behaves in that case. - Create a model in which the sections are not equally distributed and verify the correctness of the Taylor and Sachs models in this case.
- Assign a purely elastic material to one of the sections and study the behavior of the structure.
Application 1: Machine Learning Flow Rule for Hill-like plasticity¶
In this application, a machine learning algorithm is trained with data from an anisotropic Hill-like yield criterion to be used as constitutive model for anisotropic plasticity of metals.
Machine Learning flow rule¶
Machine Learning (ML) algorithms provide a great flexibility to describe aribitrary mathematical functions. At the same time they offer the possibility to handle large data sets and multi-dimensional features as input. Hence, using ML algorithms as constitutive rules for plastic material behavior offers the possibility to explicitly take into account microstructural information of the material in the constitutive modeling. Furthermore, data resulting from experiment and micromechanical simulations can be hybridized to generate training data sets. The present example is using Support Vector Classification (SVC) as yield function. The SVC algorithm is trained by using deviatoric stresses as input data and the information whether a given stress state leads to purely elastic or rather to elastic-plastic deformation of the material as result data. In this way, a ML yield function is obtained, which can determine whether a given stress state lies inside or outside of the elastic regime of the material. Furthermore, the yield locus, i.e., the hyperplane in stress space on which plastic deformation occurs, can be reconstructed from the SVC, and the gradient on this yield locus can be conveniently calculated. Therefore, the standard formulations of continuum plasticity, as the return mapping algorithm, can be applied in Finite Element Analysis (FEA) in the usual way. Thus, it is demonstrated that the new ML yield function can replace conventional FEA yield functions. In a next step, microstructural information will be considered directly as feature. Machine learning algorithms have been adopted from the scikit-learn platform (https://scikit-learn.org/stable/).
This package demonstrates the training and application of ML flow rules in FEA in form of a simple example with data synthetically produced from standard flow rules, like isotropic J2 (von Mises) and Hill-type anisotropic plasticity. It uses the pyLabFEM module. This notebook uses the matplotlib (https://matplotlib.org/) library for the visualization of results, and NumPy (http://www.numpy.org) for mathematical operations with arrays.
The scientific background of this work is desribed in mode detail in the article A. Hartmaier (2020) "Data-Oriented Constitutive Modeling of Plasticity in Metals" which is accepted for publications in Materials. The preprint is available here.
Author: Alexander Hartmaier, ICAMS, Ruhr-Universtität Bochum, Germany, March 2020
This work is licensed under a Creative
Commons Attribution-NonCommercial-ShareAlike 4.0 International License
(CC-BY-NC-SA)
The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)
Theoretical background¶
The yield function of a material is defined as \begin{equation} f = \sigma_{eq} - \sigma_y\,, \end{equation} where plastic deformation sets in at $f=0$, i.e. when the the equivalent stress $\sigma_{eq}$ equals the yield strength $\sigma_y$ of the material. The equivalent stress used here is based on the pricipal stresses $\sigma_i$ with $(i=1, 2, 3)$, as
\begin{equation} \sigma^\mathrm{J2}_{eq} = \sqrt{ \frac{1}{2}\left[ H_1 \left(\sigma_1-\sigma_2\right)^2 + H_2 \left(\sigma_2-\sigma_3\right)^2 + H_3 \left(\sigma_3-\sigma_1\right)^2 \right] } . \end{equation}In this yield function, the anisotropy of the material's flow behavior is described in a Hill-like approach for orthotropic materials. If the axes of the pricipal stresses do not match the symmetry axes of the material, the material axes and with it the parameters $H_1, H_2$ and $H_3$ must be rotated into the coordinate system of the eigenvectors of the stress tensor.
The gradient of the yield function with respect to the pricipal stresses is needed for calculating the plastic strain increments in the return mapping algorithm of continuum plasticity, and can be evaluated analytically as \begin{equation} \frac{\partial f}{\partial \sigma_1} = \frac{\partial \sigma_{eq}}{\partial \sigma_1} = \frac{\left( H_1+H_3 \right) \sigma_1 - H_1 \sigma_2 - H_3 \sigma_3}{\sigma_{eq}} \\ \frac{\partial f}{\partial \sigma_2} = \frac{\partial \sigma_{eq}}{\partial \sigma_2} = \frac{\left( H_2+H_1 \right) \sigma_2 - H_1 \sigma_1 - H_2 \sigma_3}{\sigma_{eq}} \\ \frac{\partial f}{\partial \sigma_3} = \frac{\partial \sigma_{eq}}{\partial \sigma_3} = \frac{\left( H_3+H_2 \right) \sigma_3 - H_3 \sigma_1 - H_2 \sigma_2}{\sigma_{eq}} \\ \end{equation}
Note that in the case of isotropic plasticity, i.e. $H_1=H_2=H_3=1$, the gradient takes the simple form \begin{equation} \frac{\partial f}{\partial \sigma_i} = 3 \frac{\sigma_i - p}{\sigma_{eq}} \hspace{2em} (i=1,2,3) , \end{equation} where $p = 1/3 \mbox{Tr}(\sigma)$.
Deviatoric stress space¶
Since plastic deformation in most metals does not depend on hydrostatic stress components, it is useful to transform the principal stresses from the representation as Cartesian (3D) vector $\sigma=(\sigma_1, \sigma_2, \sigma_3)$ in the principal stress space into a vector $s=(\sigma_{eq}, \theta, p)$ in the cylindrical coordinate system, with the equivalent stress $\sigma_{eq}$ representing the norm of the deviator of $\sigma$ and the polar angle $\theta$ lying in the deviatoric stress plane normal to the hydrostatic axis $p$. This coordinate transformation improves the efficiency of the training, because only 2D data for the equivalent stress and the polar angle $\theta$ need to be used as training features, whereas the hydrostatic component is disregarded. The coordinate transformation is performed by introducing a complex-valued deviatoric stress \begin{equation} \sigma'_c = \pmb{\sigma}\cdot \pmb{a} + i\,\pmb{\sigma}\cdot \pmb{b} = \sqrt{2/3}\sigma_{eq}\; e^{i\theta}, \end{equation} where $i$ is the imaginary unit, such that the polar angle is obtained as \begin{equation} \theta = \mathrm{arg}\, \sigma'_c = -i \, \ln \frac{\pmb{\sigma}\cdot \pmb{a} + i\,\pmb{\sigma}\cdot \pmb{b}} {\sqrt{2/3}\sigma_{eq}} \,, \end{equation}
where $a=(2,-1,-1)/\sqrt{6}$ (real axis) and $b=(0,1,-1)/\sqrt{2}$ (imaginary axis) are the unit vectors that span the deviatoric stress plane normal to the hydrostatic axis $c=(1,1,1)/\sqrt{3}$.
An advantage of this coordinate transformation is that the gradient of the yield function w.r.t. the cylindrical coordinates has only one non-constant component. The complete gradient w.r.t. the cylindrical coordinates reads \begin{eqnarray*} \frac{\partial f}{\partial \sigma_{eq}} &=& 1 \\ \frac{\partial f}{\partial \theta} &=& \frac{\partial \sigma_y}{\partial \theta} \\ \frac{\partial f}{\partial p} &=& 0 . \end{eqnarray*}
To transform the gradient of the flow rule from this cylindrical coordinate system back to the principle stress space, in which form it is used later on to calculate the direction of the plastic strain increments in the return mapping algorithm of the plasticity model, we introduce the Jacobian matrix for this coordinate transformation as \begin{equation} J=\frac{\partial s}{\partial \sigma} = \left( \begin{array}{ccc} \frac{\partial \sigma_{eq}}{\partial \sigma_1} & \frac{\partial \theta}{\partial \sigma_1} & \frac{\partial p}{\partial \sigma_1} \\ \frac{\partial \sigma_{eq}}{\partial \sigma_2} & \frac{\partial \theta}{\partial \sigma_2} & \frac{\partial p}{\partial \sigma_2} \\ \frac{\partial \sigma_{eq}}{\partial \sigma_3} & \frac{\partial \theta}{\partial \sigma_3} & \frac{\partial p}{\partial \sigma_3} \\ \end{array} \right) \end{equation} with ${\partial \sigma_{eq}}/{\partial \sigma_i}$ as given above, ${\partial p}/{\partial \sigma_i}=1/3$, and \begin{equation} \frac{\partial \theta}{\partial \sigma_j} = -i\, \left( \frac{\pmb{a} + i\,\pmb{b}}{\pmb{\sigma}\cdot \pmb{a} + i\,\pmb{\sigma}\cdot \pmb{b}} - \frac{3\, \pmb{\sigma}'}{\sigma_{eq}^2} \right) \hspace{2em} (j=1,2,3) . \end{equation} Finally, the gradient in the 3D principle stress space is obtained as \begin{equation} \frac{\partial f}{\partial \sigma_j} = \sum\limits_{k=1}^{3} J_{jk} \frac{\partial f}{\partial s_k} \hspace{2em} (j=1,2,3) . \end{equation}
Material definition¶
In the following, the Python class Material is invoked to be used as a material card in FEA, demonstrating the application of standard flow rules and machine learning (ML) flow rules. This class Material contains the attributes and methods to defining the elastic and plastic material behavior and to evaluate the materials constitutive behavior. Furthermore, all necessary subroutines for plotting the results are defined.
Two identical materials are defined,mat_h will be used to apply the standard Hill-like flow rules in FEA and to generate synthetical data that is used to train a ML flow rule in mat_ml. Consequently,mat_h and mat_ml should reveal identical properties, which is verified in simple FEA demonstrations for 2D plane stress cases. As reference, a further material with the same yield strength, but isotropic J2 plasticity is defined as mat_iso. Here, we consider ideal plasticity, i.e. no work hardening, and no dependence on the hydrostatic stress.
import pylabfea.model as FE
from pylabfea.material import Material
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from sklearn.metrics import r2_score
'define two elastic-plastic materials with identical yield strength and elastic properties'
E=200.e3
nu=0.3
sy = 150.
'anistropic Hill-material as reference'
mat_h = Material(name='anisotropic Hill')
mat_h.elasticity(E=E, nu=nu)
mat_h.plasticity(sy=sy, hill=[0.7,1.,1.4], drucker=0., khard=0.)
'isotropic material for ML flow rule'
mat_ml = Material(name='ML flow rule')
mat_ml.elasticity(E=E, nu=nu)
mat_ml.plasticity(sy=sy, hill=[1.,1.,1.], drucker=0., khard=0.)
print('Yield loci of anisotropic reference material and isotropic material')
ax = mat_h.plot_yield_locus(xstart=-1.5, xend=1.5, iso=True) #file='yield-locus-ref'
Create reference stress strain curves¶
The mechanical behavior of the reference material in chracterized by FEA, invoking the pyLabFEM module. Four load cases are simulated: (i) uniaxial stress in x-direction, (ii) uniaxial stress in y-direction, (iii) equibiaxial tensile strain under plane stress, and (iv) pure x-y-shear strain.
'Calculate and plot stress strain curves of reference material under various load cases'
mat_h.calc_properties(verb=False, eps=0.005, sigeps=True)
mat_h.plot_stress_strain(Hill=True) #file='sigHill-eps-ref'
Setup ML constitutive models¶
Here, the reference material mat_h with Hill-type anisotropy is used to produce training and test data for the machine learning algorithm. The yield function is represented as a step function such that Support Vector Classification (SVC) is used. In the FEA all stresses are considered in form of principal stress vectors, such that a coordinate transformation into the cylindrical coordinate system with equivalent stress and polar angle $\theta$ takes place within the training and evaluation procedure. Furthermore, to make the training process efficient, stress data for training and testing are already produced as purely deviatoric stress components, following the cylidrical stress vector definition given above.
As one cannot assume the Hill-coefficients as know for the ML material, the standard J2 equivalent stresses for isotropic materials must be used in the analysis. In consequence, the yield stress must depend on the angle $\theta$ in order to discribe the anisotropy on the flow behavior. Thus, the yield function takes the form \begin{equation} f = \sigma_{eq} - \sigma_y(\theta), \end{equation} in which the dependencies on isotropic J2 equivalent stress $\sigma_{eq}$ and angle $\theta$ are separated.
Train ML yield function¶
In the following code segment, training and test data are generated and applied for training of the SVC, and the quality test of the result. For the data and the trained ML yield function are plotted in the cylindrical stress space.
'Training and testing data for ML yield function, based on reference Material mat_h'
ndata = 36
ntest = np.maximum(20, int(ndata/10))
x_train = FE.sp_cart(mat_h.create_sig_data(ndata, extend=True)[0])
y_train = np.sign(mat_h.calc_yf(x_train))
x_test = FE.sp_cart(mat_h.create_sig_data(ntest, rand=True)[0])
y_test = np.sign(mat_h.calc_yf(x_test))
print('Plot theta vs. Hill eqiv. stress curves for reference material with known anisotropic coefficients')
ind1 = np.nonzero(y_test<0.)
ind2 = np.nonzero(y_test>=0.)
sc = FE.s_cyl(x_test, mat_h) # convert princ. stresses into cylidrical coordinates
plt.polar(sc[ind2,1],sc[ind2,0]/mat_h.sy,'.r')
plt.polar(sc[ind1,1],sc[ind1,0]/mat_h.sy,'.b')
plt.polar(np.linspace(0.,2*np.pi,36), np.ones(36), '-k', linewidth=2)
plt.legend(['test data above yield','test data below yield','Hill yield locus'], loc=(0.78,0.93))
#plt.savefig('polar-Hill.pdf', format='pdf', dpi=300)
plt.show()
'initialize and train SVC as ML yield function'
'implement ML flow rule into mat_ml'
train_sc, test_sc = mat_ml.setup_yf_SVM(x_train, y_train, x_test=x_test, y_test=y_test,
C=10, gamma=4., fs=0.3, plot=False)
y_pred_clf = mat_ml.calc_yf(x_test,pred=True)
r2_score_svm = r2_score(y_test, y_pred_clf)
print('\n-------------------------\n')
print('SVM classification fitted')
print('-------------------------\n')
print(mat_ml.svm_yf)
print("Training data points (only polar angle):", ndata,", Test data points:", ntest)
print("Training set score: {} %".format(train_sc))
print("Test set score: {} %".format(test_sc))
print("r^2 on test data : %f" % r2_score_svm)
print("Number of support vectors generated: ",len(mat_ml.svm_yf.support_vectors_))
print('Plot theta vs. J2 eqiv. stress curves for ML material')
ind1 = np.nonzero(y_test<0.)
ind2 = np.nonzero(y_test>=0.)
sc = FE.s_cyl(x_test) # convert princ. stresses into cylindrical coordinates
plt.polar(sc[ind2,1],sc[ind2,0]/mat_ml.sy,'.r')
plt.polar(sc[ind1,1],sc[ind1,0]/mat_ml.sy,'.b')
'find norm of princ. stess vector lying on yield surface'
theta = np.linspace(0.,2*np.pi,36)
snorm = FE.sp_cart(np.array([mat_ml.sy*np.ones(36)*np.sqrt(3/2), theta]).T)
x1 = fsolve(mat_ml.find_yloc, np.ones(36), args=snorm, xtol=1.e-5)
sig = snorm*np.array([x1,x1,x1]).T
s_yld = mat_ml.calc_seq(sig)
plt.polar(theta, s_yld/mat_ml.sy, '-k', linewidth=2)
plt.legend(['test data above yield','test data below yield','ML yield locus'], loc=(0.78,0.93))
#plt.savefig('polar-J2.pdf', format='pdf', dpi=300)
plt.show()
print('Plot of yield locus and training data in slices of 3D principle stress space')
mat_ml.plot_yield_locus(field=True, data=x_train, ref_mat=mat_h, trange=3.e-2,
axis1=[0,1,2], axis2=[1,2,0]) #Nmesh=300,file='yield-fct-slices'
print('Plot of trained SVM classification with test data in 2D cylindrical stress space')
xx, yy = np.meshgrid(np.linspace(-1, 1, 50),np.linspace(-1, 1, 50))
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10,8))
feat = np.c_[yy.ravel(),xx.ravel()]
#Z = mat_ml.svm_yf.predict(feat)
Z = mat_ml.svm_yf.decision_function(feat)
hl = mat_ml.plot_data(Z, ax, xx*np.pi, (yy+1.)*mat_ml.sy, c='black')
sc = FE.s_cyl(x_test)
h1 = ax.scatter(sc[:,1], sc[:,0], s=20, c=y_test, cmap=plt.cm.Paired, edgecolors='k')
ax.set_title('SVC yield function')
ax.set_xlabel(r'$\theta$ (rad)', fontsize=20)
ax.set_ylabel(r'$\sigma_{eq}$ (MPa)', fontsize=20)
ax.tick_params(axis="x", labelsize=16)
ax.tick_params(axis="y", labelsize=16)
#fig.savefig('SVM-yield-fct.pdf', format='pdf', dpi=300)
Apply trained ML flow rule in FEA¶
After the successful training of the ML yield function, it can be applied directly in FEA. Its gradients of can be calculated directly from the coefficients and the support vectors resulting from the training process. The following code segment, demonstrates the application of the ML yield function in the pyLabFEA module. First, a model with two sections is generated, where section 1 is assigned to the reference material mat_h and section 2 is assigned to mat_ml with the ML yield function.
fem=FE.Model(dim=2,planestress=True)
fem.geom([2, 2], LY=2.) # define sections in absolute lengths
print('==========================================')
print('=== FEA with reference and ML material ===')
print('==========================================')
fem.assign([mat_h, mat_ml]) # create a model with trained ML-flow rule (mat_h) and reference material (mat3)
fem.bcleft(0.)
fem.bcbot(0.)
fem.bcright(0., 'force') # free right edge
fem.bctop(0.004*fem.leny, 'disp') # apply displacement at top nodes (uniax y-stress)
fem.mesh(NX=6, NY=3)
fem.solve()
fem.calc_global()
print('Global strain: ', np.round(fem.glob['eps'],decimals=4))
print('Element strain (ref): ', np.round(fem.element[0].eps,decimals=4))
print('Element strain (ML): ', np.round(fem.element[10].eps,decimals=4))
print('Global stress: ', np.round(fem.glob['sig'],decimals=4))
print('Element stress (ref): ', np.round(fem.element[0].sig,decimals=4))
print('Element stress (ML): ', np.round(fem.element[10].sig,decimals=4))
print('Global plastic strain: ', np.round(fem.glob['epl'],decimals=5))
print('Plastic strain (ref): ', np.round(fem.element[0].epl,decimals=5))
print('Plastic strain (ML): ', np.round(fem.element[10].epl,decimals=5))
print('----------------------------')
print('Material 1 (left block): ',mat_h.msg)
print('Material 2 (right block): ',mat_ml.msg)
fem.plot('stress2',mag=1)
fem.plot('seq',mag=1)
print('Note: Different definitions of SEQ for Hill and J2.')
fem.plot('peeq',mag=1)
Note that the element stresses in both materials are very similar (top figure). However, the different definitions of the Hill equivalent stress and the J2 equivalent stress produce a numerical difference in the plot of the equivalent stresses (middle figure). The resulting plastic strains in both sections are quite similar (bottom figure).
The second example applies the same four load cases under which the reference material has been characterized and compares the results.
print('\n\n====== Stress-Strain-Curves ======')
mat_ml.calc_properties(verb=False, eps=0.005, sigeps=True, min_step=5)
mat_ml.plot_stress_strain() #file='sigJ2-eps-ML')
for sel in mat_h.prop:
hh = mat_ml.propJ2[sel]['ys']/mat_h.propJ2[sel]['ys'] - 1.
print('*** load case: ',mat_h.prop[sel]['name'])
print(' Rel. error in yield strength : ', np.round(100*hh, decimals=3),'%')
hh = np.amax(mat_ml.propJ2[sel]['peeq'])/np.amax(mat_h.propJ2[sel]['peeq']) - 1.
print('Rel. error in plastic strain at max. load : ', np.round(100*hh,decimals=3),'%')
'plot yield locus'
ax = mat_ml.plot_yield_locus(xstart=-1.5, xend=1.5, ref_mat=mat_h, field=False, Nmesh=400, fontsize=26)
print('\nPlot evolution of stresses during plastic deformation for both materials')
print('Hill-material: blue color')
print('ML-material: yellow color')
s=80
'Hill material'
stx = mat_h.sigeps['stx']['sig'][:,0:2]/mat_h.sy
sty = mat_h.sigeps['sty']['sig'][:,0:2]/mat_h.sy
et2 = mat_h.sigeps['et2']['sig'][:,0:2]/mat_h.sy
ect = mat_h.sigeps['ect']['sig'][:,0:2]/mat_h.sy
ax.scatter(stx[1:,0],stx[1:,1],s=s*2.5, c='#0000ff', edgecolors='k')
ax.scatter(sty[1:,0],sty[1:,1],s=s*2.5, c='#0000ff', edgecolors='k')
ax.scatter(et2[1:,0],et2[1:,1],s=s, c='#0000ff', edgecolors='k')
ax.scatter(ect[1:,0],ect[1:,1],s=s, c='#0000ff', edgecolors='k')
'ML material'
stx = mat_ml.sigeps['stx']['sig'][:,0:2]/mat_h.sy
sty = mat_ml.sigeps['sty']['sig'][:,0:2]/mat_h.sy
et2 = mat_ml.sigeps['et2']['sig'][:,0:2]/mat_h.sy
ect = mat_ml.sigeps['ect']['sig'][:,0:2]/mat_h.sy
ax.scatter(stx[1:,0],stx[1:,1],s=s, c='#f0ff00', edgecolors='k')
ax.scatter(sty[1:,0],sty[1:,1],s=s, c='#f0ff00', edgecolors='k')
ax.scatter(et2[1:,0],et2[1:,1],s=s, c='#f0ff00', edgecolors='k')
ax.scatter(ect[1:,0],ect[1:,1],s=s, c='#f0ff00', edgecolors='k')
#plt.savefig('comp-yield-loci.pdf', format='pdf', dpi=300)
Train ML yield function with results of 4 load cases¶
In the following, the training data is created only from the four load cases of the reference material. For all load cases $\sigma_3=0$, yet due to the data-oriented formulation of the yield function on the deviatoric plane, a full yield function can be achieved. A new material mat3 is created on which the ML training is performed.
'get principle yield stresses of 4 load cases from reference material'
sig=np.zeros((4,3))
i=0
for sel in mat_h.sigeps:
peeq = FE.eps_eq(mat_h.sigeps[sel]['epl'])
iys = np.nonzero(peeq<1.e-6) # take stress at last index of elastic regime
ys = mat_h.sigeps[sel]['sig']
sp_yld = FE.Stress(ys[iys[0][-1]]).p
sig[i,:] = sp_yld
i += 1
'mirror data in theta space: tension-compression symmetry'
sc1 = FE.s_cyl(sig)
sc2 = np.zeros((4,3))
sc2[:,0] = sc1[:,0]
sc2[:,1] = sc1[:,1]-np.pi
sc2[:,2] = sc1[:,2]
sc = np.append(sc1,sc2,axis=0)
sig = FE.sp_cart(sc)
'expand stresses into elastic and plastic regimes'
offs = 0.01
x = offs*sig
N = 23
for i in range(N):
hh = offs + (1.4-offs)*(i+1)/N
x = np.append(x, hh*sig, axis=0)
# add training points in plastic regime to avoid fallback of SVC decision fct. to zero
x = np.append(x, 2.*sig, axis=0)
x = np.append(x, 3.*sig, axis=0)
x = np.append(x, 4.*sig, axis=0)
x = np.append(x, 5.*sig, axis=0)
'result data for training of ML yield function'
y = mat_ml.calc_yf(x, pred=True)
print('Plot theta vs. J2 eqiv. stress curves for ML material')
ind1 = np.nonzero(y<0.)
ind2 = np.nonzero(y>=0.)
sc = FE.s_cyl(x) # convert princ. stresses into cylindrical coordinates
plt.polar(sc[ind2,1],sc[ind2,0]/mat_ml.sy,'.r')
plt.polar(sc[ind1,1],sc[ind1,0]/mat_ml.sy,'.b')
'find norm of princ. stess vector lying on yield surface'
theta = np.linspace(0.,2*np.pi,36)
snorm = FE.sp_cart(np.array([mat_ml.sy*np.ones(36)*np.sqrt(3/2), theta]).T)
x1 = fsolve(mat_ml.find_yloc, np.ones(36), args=snorm, xtol=1.e-5)
hs = snorm*np.array([x1,x1,x1]).T
s_yld = mat_ml.calc_seq(hs)
plt.polar(theta, s_yld/mat_ml.sy, '-k', linewidth=2)
plt.legend(['test data above yield','test data below yield','ML yield locus'], loc=(0.78,0.93))
#plt.savefig('polar-J2.pdf', format='pdf', dpi=300)
plt.show()
'define material for training with new data'
mat3 = Material(name='ML-triax')
mat3.elasticity(E=E, nu=nu)
mat3.plasticity(sy=mat_h.sy, hill=[1., 1., 1.], drucker=0., khard=0.)
'initialize and train SVC as yield function'
train_sc, test_sc = mat3.setup_yf_SVM(x, y, x_test=x_test, y_test=y_test,
C=10, gamma=4., fs=0.3, plot=False)
y_pred_clf = mat3.calc_yf(x_test,pred=True)
r2_score_svm = r2_score(y_test, y_pred_clf)
print('\n-------------------------\n')
print('SVM classification fitted')
print('-------------------------\n')
print(mat3.svm_yf)
print("Training data points (only polar angle):", ndata,", Test data points:", ntest)
print("Training set score: {} %".format(train_sc))
print("Test set score: {} %".format(test_sc))
print("r^2 on test data : %f" % r2_score_svm)
print("Number of support vectors generated: ",len(mat3.svm_yf.support_vectors_))
print('Plot of yield locus and training data in slices of 3D principle stress space')
mat3.plot_yield_locus(ref_mat=mat_h, data=x, field=True,
trange=3.e-2, axis1=[0,1,2], axis2=[1,2,0]) #Nmesh=300,file='yield-fct-slices'
Comparison of different training methods¶
The following code segment shows the same set of FE simulations as in the case of the yield function fitted to regular training data
fem=FE.Model(dim=2,planestress=True)
fem.geom([2, 2], LY=2.) # define sections in absolute lengths
print('==========================================')
print('=== FEA with reference and ML material ===')
print('==========================================')
fem.assign([mat_h, mat3]) # create a model with trained ML-flow rule (mat_h) and reference material (mat3)
fem.bcleft(0.)
fem.bcbot(0.)
fem.bcright(0., 'force') # free right edge
fem.bctop(0.004*fem.leny, 'disp') # apply displacement at top nodes (uniax y-stress)
fem.mesh(NX=6, NY=3)
fem.solve()
fem.calc_global()
print('Global strain: ', np.round(fem.glob['eps'],decimals=4))
print('Element strain (ref): ', np.round(fem.element[0].eps,decimals=4))
print('Element strain (ML): ', np.round(fem.element[10].eps,decimals=4))
print('Global stress: ', np.round(fem.glob['sig'],decimals=4))
print('Element stress (ref): ', np.round(fem.element[0].sig,decimals=4))
print('Element stress (ML): ', np.round(fem.element[10].sig,decimals=4))
print('Global plastic strain: ', np.round(fem.glob['epl'],decimals=5))
print('Plastic strain (ref): ', np.round(fem.element[0].epl,decimals=5))
print('Plastic strain (ML): ', np.round(fem.element[10].epl,decimals=5))
print('----------------------------')
print('Material 1 (left section): ',mat_h.msg)
print('Material 3 (right section): ',mat3.msg)
fem.plot('stress2',mag=1)
fem.plot('seq',mag=1)
print('Note: Different definitions of SEQ for Hill and J2.')
fem.plot('peeq',mag=1)
print('\n\n====== Stress-Strain-Curves ======')
mat3.calc_properties(verb=False, eps=0.005)
mat3.plot_stress_strain() # file='sigJ2-eps-ML'
for sel in mat_h.prop:
hh = mat_ml.propJ2[sel]['ys']/mat_h.propJ2[sel]['ys'] - 1.
print('*** load case: ',mat_h.prop[sel]['name'])
print(' Rel. error in yield strength : ', np.round(100*hh, decimals=3),'%')
hh = np.amax(mat_ml.propJ2[sel]['peeq'])/np.amax(mat_h.propJ2[sel]['peeq']) - 1.
print('Rel. error in plastic strain at max. load : ', np.round(100*hh,decimals=3),'%')
Application 2: Machine Learning Flow Rule for Tresca plasticity¶
In this application, a machine learning algorithm is trained with data from a Tresca yield criterion to be used as constitutive model for isotropic plasticity of metals.
Machine Learning flow rule¶
Machine Learning (ML) algorithms provide a great flexibility to describe aribitrary mathematical functions. At the same time they offer the possibility to handle large data sets and multi-dimensional features as input. Hence, using ML algorithms as constitutive rules for plastic material behavior offers the possibility to explicitly take into account microstructural information of the material in the constitutive modeling. Furthermore, data resulting from experiment and micromechanical simulations can be hybridized to generate training data sets. The present example is using Support Vector Classification (SVC) as yield function. The SVC algorithm is trained by using deviatoric stresses as input data and the information whether a given stress state leads to purely elastic or rather to elastic-plastic deformation of the material as result data. In this way, a ML yield function is obtained, which can determine whether a given stress state lies inside or outside of the elastic regime of the material. Furthermore, the yield locus, i.e., the hyperplane in stress space on which plastic deformation occurs, can be reconstructed from the SVC, and the gradient on this yield locus can be conveniently calculated. Therefore, the standard formulations of continuum plasticity, as the return mapping algorithm, can be applied in Finite Element Analysis (FEA) in the usual way. Thus, it is demonstrated that the new ML yield function can replace conventional FEA yield functions. In a next step, microstructural information will be considered directly as feature. Machine learning algorithms have been adopted from the scikit-learn platform (https://scikit-learn.org/stable/).
This package demonstrates the training and application of ML flow rules in FEA in form of a simple example with data synthetically produced from standard flow rules, like isotropic J2 (von Mises) and Hill-type anisotropic plasticity. It uses the pyLabFEM module. This notebook uses the matplotlib (https://matplotlib.org/) library for the visualization of results, and NumPy (http://www.numpy.org) for mathematical operations with arrays.
The scientific background of this work is desribed in mode detail in the article A. Hartmaier (2020) "Data-Oriented Constitutive Modeling of Plasticity in Metals" which is accepted for publications in Materials. The preprint is available here.
Author: Alexander Hartmaier, ICAMS, Ruhr-Universtität Bochum, Germany, March 2020
This work is licensed under a Creative
Commons Attribution-NonCommercial-ShareAlike 4.0 International License
(CC-BY-NC-SA)
The pyLabFEA package comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under the conditions of the GNU General Public License (GPLv3)
Theoretical background¶
A more detailed description of the theoretical background can be found in the notebook introducing ML flow rules for materials with Hill-like anisotropic plasticity. The online docummentation of the entire pyLabFEA package is available here. In this notebook, only the fundamental quantities are introduced. The yield function of a material is defined as \begin{equation} f = \sigma_{eq} - \sigma_y\,, \end{equation} where plastic deformation sets in at $f=0$, i.e. when the the equivalent stress $\sigma_{eq}$ equals the yield strength $\sigma_y$ of the material. The Trseca equivalent stress used here is based on the pricipal stresses $\sigma_\mathrm{I} \ge \sigma_\mathrm{II} \ge \sigma_\mathrm{III}$, as
\begin{equation} \sigma^\mathrm{Tresca}_{eq} = \frac{\sigma_\mathrm{I} - \sigma_\mathrm{III}}{2} . \end{equation}Material definition¶
In the following, the Python class Material is invoked to be used as a material card in FEA, demonstrating the application of standard flow rules and machine learning (ML) flow rules. This class Material contains the attributes and methods to defining the elastic and plastic material behavior and to evaluate the materials constitutive behavior. Furthermore, all necessary subroutines for plotting the results are defined.
Two identical materials are defined,mat_t will be used to define the Tresca flow rules for generating synthetical data that is used to train a ML flow rule in mat_ml. As reference, a further material with the same yield strength, but isotropic J2 plasticity is defined as mat_iso. Here, we consider ideal plasticity, i.e. no work hardening, and no dependence on the hydrostatic stress.
import pylabfea.pyLabFEM as FE
from pylabfea.pyLabMaterial import Material
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from sklearn.metrics import r2_score
'define two elastic-plastic materials with identical yield strength and elastic properties'
E=200.e3
nu=0.3
sy = 150.
'Tresca-material as training model'
mat_t = Material(name='Tresca')
mat_t.elasticity(E=E, nu=nu)
mat_t.plasticity(sy=sy, Tresca=True)
mat_iso = Material(name='isotropic J2')
mat_iso.elasticity(E=E, nu=nu)
mat_iso.plasticity(sy=sy)
'isotropic material as basis for ML flow rule'
mat_ml = Material(name='ML flow rule')
mat_ml.elasticity(E=E, nu=nu)
mat_ml.plasticity(sy=sy, hill=[1.,1.,1.], drucker=0., khard=0.)
print('Yield loci of anisotropic reference material and isotropic material')
ax = mat_t.plot_yield_locus(xstart=-1.5, xend=1.5, ref_mat=mat_iso, field=True) #file='yield-locus-ref'
Setup ML constitutive models¶
Here, the reference material mat_t with the Tresca yield function is used to produce training and test data for the machine learning algorithm. The yield function is represented as a step function such that Support Vector Classification (SVC) is used.
Cylindrical coordinates in deviatoric stress space¶
Since plastic deformation in most metals does not depend on hydrostatic stress components, the training process can be simplified by transforming the principal stresses from the representation as Cartesian (3D) vector $\sigma=(\sigma_1, \sigma_2, \sigma_3)$ in the principal stress space into a vector $s=(\sigma_{eq}, \theta, p)$ in the cylindrical coordinate system, with the equivalent stress $\sigma_{eq}$ representing the norm of the deviator of $\sigma$ and the polar angle $\theta$ lying in the deviatoric stress plane normal to the hydrostatic axis $p$. This coordinate transformation improves the efficiency of the training, because only 2D data for the equivalent stress and the polar angle $\theta$ need to be used as training features, whereas the hydrostatic component is disregarded. The coordinate transformation is performed by introducing a complex-valued deviatoric stress \begin{equation} \sigma'_c = \pmb{\sigma}\cdot \pmb{a} + i\,\pmb{\sigma}\cdot \pmb{b} = \sqrt{2/3}\sigma_{eq}\; e^{i\theta}, \end{equation} where $i$ is the imaginary unit, such that the polar angle is obtained as \begin{equation} \theta = \mathrm{arg}\, \sigma'_c = -i \, \ln \frac{\pmb{\sigma}\cdot \pmb{a} + i\,\pmb{\sigma}\cdot \pmb{b}} {\sqrt{2/3}\sigma_{eq}} \,, \end{equation}
where $a=(2,-1,-1)/\sqrt{6}$ (real axis) and $b=(0,1,-1)/\sqrt{2}$ (imaginary axis) are the unit vectors that span the deviatoric stress plane normal to the hydrostatic axis $c=(1,1,1)/\sqrt{3}$.
Train ML yield function¶
In the following code segment, training and test data are generated and applied for training of the SVC, and the quality test of the result. To make the training process efficient, stress data for training and testing are already produced as purely deviatoric stress components, following the cylidrical stress vector definition given above. For the data and the trained ML yield function are plotted in the cylindrical stress space.
'Create training data in deviatoric stress space for components seq and theta'
def create_data(N, mat, extend=False, rand=False):
# create stresses along unit circle normal to hydrostatic axis
if not rand:
theta = np.linspace(-np.pi, np.pi, N)
else:
theta = 2.*(np.random.rand(N)-0.5)*np.pi
sig = np.array([np.ones(N), theta]).T
'expand stresses into elastic and plastic regimes'
'construct elastic data points'
xt = sig*np.array([0.01,1.])
offs = 0.05
nl = 15
for i in range(nl):
hh = np.array([1 - (i/nl)**1.1, 1.]) # down-scaling of equiv. stresses in elastic regime
xt = np.append(xt, hh*sig, axis=0)
'construct plastic data points'
#xt = np.append(xt, sig, axis=0)
offs = 1.01
for i in range(nl):
hh = np.array([offs + (i/nl)**1.1, 1.])
xt = np.append(xt, hh*sig, axis=0)
'add further training points far in plastic regime to avoid fallback of SVC decision fct. to zero'
if extend:
for i in range(9):
hh = np.array([i*0.5+2., 1.])
xt = np.append(xt, hh*sig, axis=0)
'result data for ML yield function (only sign is considered)'
x = FE.sp_cart(xt)
y = np.sign(mat.calc_yf(x*mat.sy))
return x,y
'Training and testing data for ML yield function, based on reference Material mat_t'
ndata = 600
ntest = np.maximum(20, int(ndata/10))
x_train, y_train = create_data(ndata, mat_t, extend=True)
x_test, y_test = create_data(ntest, mat_t, rand=False)
x_train *= mat_t.sy
x_test *= mat_t.sy
print('Plot theta vs. Tresca eqiv. stress curves for reference material')
ind_e = np.nonzero(y_test<0.)[0]
ind_pl = np.nonzero(y_test>=0.)[0]
sc = FE.s_cyl(x_test, mat_t) # convert princ. stresses into cylidrical coordinates
plt.polar(sc[ind_pl,1],sc[ind_pl,0]/mat_t.sy,'.r')
plt.polar(sc[ind_e,1],sc[ind_e,0] /mat_t.sy,'.b')
'find norm of princ. stess vector lying on yield surface'
theta = np.linspace(-np.pi,np.pi,36)
snorm = FE.sp_cart(np.array([mat_t.sy*np.ones(36)*np.sqrt(3/2), theta]).T)
x1 = fsolve(mat_t.find_yloc, np.ones(36), args=snorm, xtol=1.e-5)
sig = snorm*np.array([x1,x1,x1]).T
s_yld = mat_t.calc_seq(sig)
plt.polar(theta, s_yld/mat_t.sy, '-k', linewidth=2)
plt.legend(['training data (plastic)','training data (elastic)','Tresca yield locus'], loc=(0.78,0.93))
#plt.savefig('polar-Hill.pdf', format='pdf', dpi=300)
plt.show()
'initialize and train SVC as ML yield function'
'implement ML flow rule into mat_ml'
train_sc, test_sc = mat_ml.setup_yf_SVM(x_train, y_train, x_test=x_test, y_test=y_test,
C=50, gamma=9., fs=0.4, plot=False)
y_pred_clf = mat_ml.calc_yf(x_test,pred=True)
r2_score_svm = r2_score(y_test, y_pred_clf)
print('\n-------------------------\n')
print('SVM classification fitted')
print('-------------------------\n')
print(mat_ml.svm_yf)
print("Training data points (only polar angle):", ndata,", Test data points:", ntest)
print("Training set score: {} %".format(train_sc))
print("Test set score: {} %".format(test_sc))
print("r^2 on test data : %f" % r2_score_svm)
print("Number of support vectors generated: ",len(mat_ml.svm_yf.support_vectors_))
print('Plot theta vs. J2 eqiv. stress curves for ML material')
ind1 = np.nonzero(y_test<0.)
ind2 = np.nonzero(y_test>=0.)
sc = FE.s_cyl(x_test) # convert princ. stresses into cylindrical coordinates
plt.polar(sc[ind2,1],sc[ind2,0]/mat_ml.sy,'.r')
plt.polar(sc[ind1,1],sc[ind1,0]/mat_ml.sy,'.b')
'find norm of princ. stess vector lying on yield surface'
theta = np.linspace(-np.pi,np.pi,36)
snorm = FE.sp_cart(np.array([mat_ml.sy*np.ones(36)*np.sqrt(3/2), theta]).T)
x1 = fsolve(mat_ml.find_yloc, 0.8*np.ones(36), args=snorm, xtol=1.e-5)
sig = snorm*np.array([x1,x1,x1]).T
s_yld = mat_ml.calc_seq(sig)
plt.polar(theta, s_yld/mat_ml.sy, '-k', linewidth=2)
plt.legend(['test data above yield','test data below yield','ML yield locus'], loc=(0.78,0.93))
#plt.savefig('polar-J2.pdf', format='pdf', dpi=300)
plt.show()
print('Plot of yield locus and training data in slices of 3D principle stress space')
mat_ml.plot_yield_locus(field=True, data=x_train, ref_mat=mat_t, trange=3.e-2,
axis1=[0,1,2], axis2=[1,2,0]) #, Nmesh=400, file='Tresca-yfct-slices')
print('Plot of trained SVM classification with test data in 2D cylindrical stress space')
xx, yy = np.meshgrid(np.linspace(-1, 1, 50),np.linspace(-1, 1, 50))
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10,8))
feat = np.c_[yy.ravel(),xx.ravel()]
#Z = mat_ml.svm_yf.predict(feat)
Z = mat_ml.svm_yf.decision_function(feat)
hl = mat_ml.plot_data(Z, ax, xx*np.pi, (yy+1.)*mat_ml.sy, c='black')
sc = FE.s_cyl(x_test)
h1 = ax.scatter(sc[:,1], sc[:,0], s=20, c=y_test, cmap=plt.cm.Paired, edgecolors='k')
#ax.set_title('SVC yield function')
ax.set_xlabel(r'$\theta$ (rad)', fontsize=20)
ax.set_ylabel(r'$\sigma_{eq}$ (MPa)', fontsize=20)
ax.tick_params(axis="x", labelsize=16)
ax.tick_params(axis="y", labelsize=16)
#fig.savefig('SVM-yield-fct_Tresca.pdf', format='pdf', dpi=300)
Apply trained ML flow rule in FEA¶
After the successful training of the ML yield function, it can be applied directly in FEA. Its gradients of can be calculated directly from the coefficients and the support vectors resulting from the training process. The following code segment, demonstrates the application of the ML yield function in the pyLabFEA module. First, a model with two sections is generates, where section 1 is assigned to the reference material mat_t and section 2 is assigned to mat_ml with the ML yield function.
fem=FE.Model(dim=2,planestress=True)
fem.geom([2, 2], LY=2.) # define sections in absolute lengths
print('==========================================')
print('=== FEA with reference and ML material ===')
print('==========================================')
fem.assign([mat_iso, mat_ml]) # create a model with trained ML-flow rule (mat_h) and reference material (mat3)
fem.bcleft(0.)
fem.bcbot(0.)
fem.bctop(0.004*fem.leny, 'disp') # apply displacement at top nodes (uniax y-stress)
fem.mesh(NX=6, NY=3)
fem.solve()
fem.calc_global()
print('Global strain: ', np.round(fem.glob['eps'],decimals=4))
print('Element strain (ref): ', np.round(fem.element[0].eps,decimals=4))
print('Element strain (ML): ', np.round(fem.element[10].eps,decimals=4))
print('Global stress: ', np.round(fem.glob['sig'],decimals=4))
print('Element stress (ref): ', np.round(fem.element[0].sig,decimals=4))
print('Element stress (ML): ', np.round(fem.element[10].sig,decimals=4))
print('Global plastic strain: ', np.round(fem.glob['epl'],decimals=5))
print('Plastic strain (ref): ', np.round(fem.element[0].epl,decimals=5))
print('Plastic strain (ML): ', np.round(fem.element[10].epl,decimals=5))
print('----------------------------')
print('Material 1 (left block): ',mat_iso.msg)
print('Material 2 (right block): ',mat_ml.msg)
fem.plot('seq',mag=1)
fem.plot('peeq',mag=1)
Note that the J2 equivalent stresses in both materials are somewhat different althuogh under a uniaxial stress, the Tresca and J2 yield stresse should be the same. The difference results from the rounding-off of the sharp corners of the Tresca yield function by the ML yield function. The resulting plastic strains in both sections are still quite similar (bottom figure).
The second example applies the same four load cases under which the reference material has been characterized and compares the results.
print('\n\n====== Stress-Strain-Curves ======')
mat_ml.calc_properties(verb=False, eps=0.005, sigeps=True)
mat_ml.plot_stress_strain()
print('Yield loci and flow stresses in sig_1-sig_2 principle stress space')
ax = mat_ml.plot_yield_locus(xstart=-1.5, xend=1.5, ref_mat=mat_t, iso=True, field=False, Nmesh=400, fontsize=26)
s=100
'ML material'
stx = mat_ml.sigeps['stx']['sig'][:,0:2]/mat_ml.sy
sty = mat_ml.sigeps['sty']['sig'][:,0:2]/mat_ml.sy
et2 = mat_ml.sigeps['et2']['sig'][:,0:2]/mat_ml.sy
ect = mat_ml.sigeps['ect']['sig'][:,0:2]/mat_ml.sy
ax.scatter(stx[1:,0],stx[1:,1],s=s, c='#f0ff00', edgecolors='k')
ax.scatter(sty[1:,0],sty[1:,1],s=s, c='#f0ff00', edgecolors='k')
ax.scatter(et2[1:,0],et2[1:,1],s=s, c='#f0ff00', edgecolors='k')
ax.scatter(ect[1:,0],ect[1:,1],s=s, c='#f0ff00', edgecolors='k')
plt.show()
sect = mat_t.calc_seq(mat_ml.sigeps['ect']['sig'][1,0:3])
set2 = mat_t.calc_seq(mat_ml.sigeps['et2']['sig'][1,0:3])
stx = mat_t.calc_seq(mat_ml.sigeps['stx']['sig'][1,0:3])
sty = mat_t.calc_seq(mat_ml.sigeps['sty']['sig'][1,0:3])
print('\n====== Tresca equivalent stress at yielding ======')
print('Nominal yield strength: ',mat_ml.sy,'MPa')
print('load case: ',mat_ml.prop['stx']['name'], '; Tresca equivalent stress: ',stx.round(decimals=2),'MPa')
print('load case: ',mat_ml.prop['sty']['name'], '; Tresca equivalent stress: ',sty.round(decimals=2),'MPa')
print('load case: ',mat_ml.prop['et2']['name'], '; Tresca equivalent stress: ',set2.round(decimals=2),'MPa')
print('load case: ',mat_ml.prop['ect']['name'], '; Tresca equivalent stress: ',sect.round(decimals=2),'MPa')
#plt.savefig('comp-yield-loci_Tresca.pdf', format='pdf', dpi=300)
It is seen that the flow stresses of the ML yield function stay on the yield locus, as it is expect for ideal plasticity.